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I.  INTRODUCTION 

This  report  is  the  Systems,  Science  and  Software  (S3) 
contribution  to  a joint  research  project  carried  out  in 
collaboration  with  Applied  Theory,  Incorporated  (ATI)  and 
Pacific  Sierra  Research  (PSR) . A series  of  cratering  calcula- 
tions were  done  by  ATI  and  the  results  were  transmitted  to 
S3  for  processing.  The  S3  contribution  has  been  to  compute 
theoretical  body  and  surface  wave  records  from  the  ATI  data 
and  from  these  to  obtain  m^  and  Mg.  The  results  are  to  be 
analyzed  by  PSR. 

The  computational  procedure  employed  at  S3  is  to 
first  obtain  an  equivalent  elastic  source  representation  of 
the  cratering  explosions.  Having  an  elastic  source  repre- 
sentation, elastic  wave  propagation  methods  can  be  employed 
to  compute  theoretical  seismograms.  The  key  step  in  this 
procedure  is  clearly  the  interpretation  of  the  raw  data  in 
terms  of  an  equivalent  elastic  source. 

The  computation  of  an  equivalent  elastic  source  from 
the  output  of  finite  difference  calculations  has  been  ac- 
complished at  S3  for  a number  of  complex  explosion  and  earth- 
quake sources.  However,  all  of  our  experience  has  been  with 
sources  that  were  entirely  surrounded  by  an  elastic  material. 
The  cratering  sources  then  provided  a special  problem.  In 
fact,  our  technique  is  not  rigorously  valid  for  half-space 
sources  of  this  kind.  Therefore,  it  is  quite  important  to 
understand  the  approximations  that  must  be  employed  and  their 
potential  effect  on  the  results. 

In  analyzing  the  approximations  made  we  conclude  that 
in  every  case  the  effects  are  much  more  severe  for  the  S 
waves  then  for  the  P waves.  In  fact,  the  short  period  P 
waves  seem  to  be  handled  very  well  by  our  technique.  Since 
the  S waves  have  almost  no  effect  on  the  body  wave  seismo- 
grams, we  have  considerable  confidence  in  the  it^  values 
generated . 
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In  order  to  compute  Mg  we  have  to  accurately  compute 
the  source  at  periods  of  about  20  seconds.  The  finite  dif- 
ference calculations  were  terminated  at  2.5  seconds.  There- 
fore, our  Mg  calculations  depend  on  the  source  havina  reached 
a static  state.  This  is  much  closer  to  being  true  for  the  P 
waves  than  for  the  S waves.  Fortunately,  we  find  that  for 
this  case  the  S waves  have  to  be  quite  a lot  larger  than  the 
P waves  to  substantially  perturb  Mg,  so  our  errors  in  the 
long  period  S waves  may  not  be  too  important.  Still,  we  have 

less  confidence  in  the  M values  than  in  the  m,  values. 

s b 

Nevertheless,  we  are  treating  all  calculations  the  same  and 

we  can  reasonably  suppose  that  trends  are  properly  reflected 

in  the  M data, 
s 

We  also  computed  the  seismic  waves  for  the  ejecta 
fallback  for  several  of  the  sources.  We  find  that  this 
contribution  to  the  body  and  surface  waves  is  too  small  to 
be  of  any  significance  for  m^  and  Mg  measurements. 

The  report  is  organized  in  nine  sections  and  five 
appendices.  In  Section  II  we  outline  the  computational 
procedures.  We  also  summarize  the  content  of  the  appendices 
where  detailed  exposition  of  some  aspects  of  the  theory  and  to 
application  is  given.  In  Section  III  we  briefly  describe  the 
fourteen  finite  difference  calculations  to  be  studied. 

Twelve  of  these  are  cratering  calculations  and  two  are  for 
spherically  symmetric  contained  explosions.  In  Section  IV 
we  discuss  the  equivalent  elastic  source  representation  of 
the  cratering  calculations.  The  approximations  made  and 
their  effect  on  the  solution  is  indicated  in  some  detail  in 
this  section.  The  far-field  displacement  spectra  are 


extracted  from  the  source  calculations  and  displayed  in 
Section  V.  Sections  VI  and  VII  give  the  m^  and  Mg  results. 
In  Section  VII  the  seismic  waves  generated  by  the  ejecta 
fallback  are  discussed.  Finally,  in  Section  IX  the  sources 


are  scaled  to  37.5  and  600  kt  and  the  m^  and  M are  computed 
for  each  of  the  fourteen  sources. 
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II.  COMPUTATIONAL  PROCEDURE 


In  this  section  we  outline  the  computational  procedure 
used  to  compute  Mg  and  m^  for  the  ATI  cratering  calculations . 
The  main  steps  in  our  procedure  may  be  listed  as  follows : 

1.  The  ATI  performed  cratering  calculations  were 
carried  into  the  regime  where  the  material  res- 
ponse is  approximately  linearly  elastic.  A tape 
was  then  prepared  containing  the  time  histories 
of  the  divergence  and  curl  (V»u(t),  Vxu(t))  of 
the  displacement  field  on  a radius,  denoted  the 
elastic  radius,  centered  at  ground  zero.  The 
geometry  and  coordinate  system  are  shown  in 
Figure  2.1. 

2.  The  divergence  and  curl  are  expanded  in  a series 
of  spherical  harmonics  to  obtain  an  equivalent 
elastic  source.  The  procedure  is  formally  that 
described  by  Bache  and  Harkrider  [1976]  for 
sources  in  a whole  space.  However,  the  presence 
of  a free  surface  requires  a number  of  assumptions 
that,  to  some  extent,  control  the  solution.  The 
extent  of  this  effect  is  discussed  in  some  detail 
in  this  report.  However,  the  most  important 
point  is  that  all  the  cratering  calculations  are 
treated  the  same  way  and  the  relative  values  of 

m^  and  Mg  should  be  preserved. 

3.  Using  the  equivalent  elastic  source,  synthetic 
seismograms  are  computed  for  body  waves  and  sur- 
face waves.  For  body  waves  the  pertinent  references 
are  Bache  and  Harkrider  [1976]  and  Bache,  et  al. 
[1976] . For  surface  waves  we  use  the  method  of 
Harkrider  which  has  been  described  in  numerous 
publications;  e.g.,  Harkrider  [1964],  with  certain 
modifications  indicated  in  Section  VII. 
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Figure  2.1.  The  geometry  and  coordinate  system  for  the  crater 
ing  calculations.  The  solution  is  independent  of 
the  azimuthal  coordinate. 
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4.  The  values  of  m^  and  Mg  are  strongly  dependent  on 
the  crust  and  upper  mantle  models  used  in  the  syn- 
thetic seismogram  calculations.  However,  the  same 
models  are  used  in  all  calculations  and  the  rela- 
tive values  should  be  insensitive  to  these  models. 

5.  The  free  surface  stresses  due  to  the  ejecta  fall- 
back were  also  provided  to  us  by  ATI.  These  data 
were  analyzed  to  determine  the  effect  on  the  far- 
field  body  and  surface  waves. 

Further  descriptions  of  the  computational  procedure 
are  contained  in  subsequent  sections  of  the  main  body  of  the 
report  where  the  equivalent  elastic  source  representation, 
the  far-field  displacement  spectra  and  the  teleseismic  wave- 
forms are  discussed.  Also  included  in  the  report  are  a 
number  of  Appendices  which  treat  various  aspects  of  our  work 
in  considerable  detail.  The  content  of  these  appendices  is 
briefly  summarized  below. 

Appendix  A:  Theory  of  the  Equivalent  Elastic  Source  Repre- 

sentation 

This  appendix  is  a reproduction  of  a section  of  a 
paper  by  Bache  and  Harkrider  [1976]  . It  outlines  the  mathe- 
matical formulation  for  representing  an  arbitrary  volume 
source  in  an  elastic,  homogeneous  and  isotropic  space  in 
terms  of  an  expansion  in  spherical  harmonics.  The  expansion 
coefficients,  the  multipole -coefficients , provide  an  equiva- 
lent elastic  source  representation. 

Appendix  B:  Application  of  the  Multipolar  Expansion  Techni- 

gue  to  Cratering  Explosions 

The  procedures  of  Appendix  A,  appropriately  modified, 
can  give  a satisfactory  representation  of  the  output  of  the 
cratering  calculations.  The  approximations  required  are 
discussed  in  Appendix  B. 
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Appendix  C:  Effect  of  Source  Symmetry  on  the  Multipolar 

Source 

The  most  important  approximation  is  in  the  choice  of 
the  symmetry  to  be  imposed  on  the  computed  radiation  field. 
The  possibilities  are  discussed  in  Appendix  C. 

Appendix  D;  Detailed  Description  of  a Typical  Calculation 
of:  the  Multipolar  Source  Representation 

We  go  through  a typical  calculation  step-by-step  and 
indicate  the  operations  carried  out. 

Appendix  E:  The  Seismic  Waves  Due  to  a Stress  Distribution 

Applied  at  the  Surface  of  a Multilayered  Half- 
space 

In  order  to  compute  m^  and  M for  the  ejecta  impact 
portion  of  the  cratering  calculations,  it  was  necessary  to 
extend  our  theoretical  results  to  include  the  case  of  a 
distribution  load  on  the  surface.  The  theory  is  developed 
in  Appendix  E. 
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III.  DESCRIPTION  OF  THE  CRATERING  CALCULATIONS 


In  order  to  carry  out  our  analyses  of  the  cratering 
calculations,  it  was  necessary  to  have  some  information 
about  the  source  depth  and  the  local  material  properties. 

This  information , together  with  an  identifier  for  each 
calculation  for  use  in  subsequent  sections,  is  summarized 
here . 

We  will  be  describing  fourteen  source  calculations  in 
three  emplacement  materials.  Two  of  these  are  one-dimen- 
sional calculations  for  contained  explosions  in  a homogeneous 
whole  space  and  provide  benchmark  cases  for  measuring  the 
effect  of  the  cratering.  The  other  twelve  are  cratering 
calculations.  The  important  parameters  characterizing  the 
calculations  are  summarized  in  Table  3.1. 

We  see  from  the  table  that  the  calculations  have  the 
potential  to  help  us  understand: 

1.  The  effect  of  burial  depth  on  the  teleseismic 
signature  of  cratering  explosions. 

2.  The  effect  of  emplacement  material  on  the  signal 
from  cratering  explosions. 

3.  The  difference  between  cratering  and  contained 
shots  in  the  same  material. 

Considering  the  approximations  made  in  our  calculations  (ex- 
clusive of  any  difficulties  with  the  finite  difference  cal- 
culations) , the  above  are  listed  in  order  of  the  amount  of 
confidence  we  place  in  the  results.  That  is,  we  are  most 
optimistic  about  our  conclusions  regarding  differences  in 
the  source  coupling  between  events  that  differ  only  by  the 
burial  depth.  On  the  other  hand,  we  are  least  confident 
in  conclusions  about  the  differences  between  cratering  and 
contained  explosions. 
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TABLE  3.1 

PARAMETERS  DESCRIBING  THE  ATI  CALCULATIONS 


Identifier 

Material 

Depth  (km) 

a(km/sec) 

6 (km/sec) 

p (km/ sec) 

1 

Granite 

0.159 

4.402 

2.54 

2.661 

2 

II 

0.207 

4.402 

2.54 

2.661 

3 

II 

0.253 

4.406 

2.542 

2.661 

4 

Dry  Sandstone 

0.159 

2.322 

1.740 

2.30 

5 

II 

0.207 

2.825 

1.743 

2.30 

6 

It 

0.253 

2.828 

1.744 

2.30 

7 

Low  Strength 

0.207 

2.836 

1.755 

2.30 

8 

Wet  Sandstone 

0.053 

2.620 

1.509 

2.40 

9 

II 

0.159 

2.624 

1.592 

2.40 

10 

II 

0.207 

2.614 

1.513 

2.40 

11 

II 

0.253 

2.613 

1.517 

2.40 

12 

II 

0.531 

2.619 

1.519 

2.40 

13 

Granite 

Spherically 

Symmetric 

4.239 

2.448 

2.661 

14 

Wet  Sandstone 

Spherically 

Symmetric 

2.619 

1.530 

2.40 
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IV.  EQUIVALENT  ELASTIC  SOURCE  REPRESENTATION 


4 . 1 MULTIPOLE  COEFFICIENTS  IN  THE  TIME  DOMAIN 

The  procedure  used  to  obtain  the  equivalent  elastic 
source  representation  for  the  cratering  calculations  is  des- 
cribed in  detail  in  the  Appendices.  The  basic  equations 
written  in  the  frequency  domain  are  given  by  A.1-A.5  of 
Appendix  A.  In  our  implementation,  we  actually  use  the  time 
domain  analogues  of  these  equations.  That  is,  the  equation 
of  motion  is  written 


= a2  V X{4)  - 2B2  7 x x 


(4.1) 


where  the  Cartesian  potentials  are  defined  by 
X4  = V • u 

X = j V x u. 


(4.2) 


The  potentials  are  expanded  in  spherical  eigenfunctions 
as  follows 


(j) 


(R, 


A^ ^ (R,  t)  cos  m<J> 

Jtm 


£=0  m=»0 
+ B ^ (R,  t)  sin  m<|> 


(4.3) 


m 


P™(cos6)  , j = 1,2, 3, 4 


Note  that  these  equations  are  identical  to  A.l,  A. 2,  A. 4,  if 


(to) 


hj21  (K  ^ R) 


/ Ata(R't)e‘t“t  dt- 

— oo 


(4.4) 
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and  similarly  for  (w) . Note  that  (4.4)  is  simply  a 

Fourier  transformation  normalized  by  the  Hankel  function. 

In  computing  the  multipolar  representation  of  the 
cratering  calculations  we  assume  axisymmetry  and  vertical 
antisymmetry,  the  case  described  in  Section  C.2.3  of 
Appendix  C.  Then  the  multipole  coefficients  in  the  time 
domain  are  actually  computed  from  C.32-C.34.  For  example, 
for  the  P-wave  portion  of  the  field  the  nonzero  terms  are 

tt/2 

(R,t)  = (24+1)  / x(4)(R,6,t)  pJ  (cosG)  sin6  d0 , 

0 

l = 1,3,5. . . , (4.5) 

(4 ) 

and  the  x (R,0,t)  is  the  time  history  of  the  divergence 
of  the  displacement  field  given  on  the  quarter-circle  of 
radius  R0  (see  Figure  2.1). 

It  turns  out  that  the  dominant  term  for  the  m^  and  Mg 

determinations  is  the  leading  term  in  the  expansion  of  the 

(4 ) 

P-wave  portion  of  the  field;  that  is,  (R,t).  The  other 

terms  provide  higher  order  corrections  as  will  be  explained 

in  later  sections.  For  now,  let  us  examine  the  "dipole-P" 

term  as  it  is  a fair  representation  of  the  equivalent  elastic 

(4 ) 

source.  In  Figures  4.1  through  4.3  the  computed  A^Q  (R,t) 
is  plotted  for  each  of  the  twelve  cratering  calculations. 

The  multipole  coefficient  is  dimensionless  and  the  amplitude 
depends  on  the  radius,  R. 

For  each  calculation  we  indicate  the  theoretical  ar- 
rival time  for  the  P,  S and  Rayleigh  (R)  waves.  For  the 
P wave  this  time  is  simply  (R-h)/a  since  the  shortest  path 
from  the  source  to  the  elastic  radius  is  straight  down.  The 
same  is  done  for  the  S wave.  The  theoretical  arrival  time 

for  the  Rayleigh  wave  at  the  free  surface  is  R/V_.»  where 

R 

the  Rayleigh  wave  velocity,  V = 0.91926.  This  is  the  time 
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gure  4.3.  The  dipole-P  multipole  coefficient  for  the  five  cratering  calculations  in 
wet  sandstone.  For  Case  8,  R = 1.3  km  and  R = 1.2  km  for  the  remaining 
four  calculations. 
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indicated  on  the  plots.  However,  the  Rayleigh  wave  arrives 
earlier  at  depth,  though  its  amplitude  decays  exponentially 
with  depth.  For  example,  at  30°  down  from  the  free  surface 
the  arrival  time  is  cos  30°  or  0.866  times  that  shown.  This 
time  turns  out  to  be  close  to  the  S wave  arrival  time  indi- 
cated so  the  Rayleigh  wave  arrival  can  be  thought  to  be 
somewhere  in  the  zone  between  the  two.  Otherwise,  no  energy 
is  expected  at  the  S wave  arrival  time  as  no  S waves  should 
be  propagated  directly  from  the  source.  Further,  the  diver- 
gence should  be  transparent  to  S wave  radiation. 

Examining  the  dipole-P  term  for  the  granite  calcula- 
tions, the  theoretical  P wave  arrival  is  an  elastic  travel 
time  and  is  later  than  the  actual  energy  arrival.  This  is 
expected  because  of  the  higher  shock  wave  velocities  in  the 
nonlinear  regime  near  the  explosion.  There  is  no  indication 
of  the  Rayleigh  wave  having  a substantial  influence  on  the 
record.  Even  if  the  Rayleigh  wave  displacements  were  large, 
we  would  expect  them  to  perturb  the  curl  more  strongly 
than  the  divergence.  Finally,  we  observe  that  after  the 
passage  of  the  main  pulse,  the  time  series  approaches  zero 
as  a static  limit  and  oscillates  about  this  value  by  small, 
and  probably  insignificant,  amounts. 

For  the  dry  sandstone  case  (Figure  4.2)  the  dipole-P 
coefficient  behaves  differently  than  for  the  granite  after 
the  first  peak  and  trough.  Once  again,  it  seems  to  settle 
to  a static  value  at  about  the  time  marked  R.  However,  this 
static  value  is  clearly  negative  and  the  oscillations  are 
considerably  larger  compared  to  the  peak  values. 

The  behavior  of  the  dipole-P  coefficient  is  least 
consistent  for  the  wet  sandstone  cases  shown  in  Figure  4.3. 
The  Case  8 term  behaves  like  those  for  granite  while  Cases 
9 and  10  look  more  like  the  dry  sandstone  coefficients.  In 
Cases  11  and  12  is  the  only  indication  of  possible  difficulty 
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associated  with  the  presence  of  Rayleigh  wave  energy.  In 
the  other  cases  the  static  limit  is  essentially  attained  be- 
fore the  Rayleigh  wave  arrival  time  but  this  is  not  true  for 
the  last  two  cases. 

Finally,  we  should  mention  that  some  of  the  sudden 
jumps  that  occur  in  these  time  series,  most  prominently  those 
in  the  first  peak  of  Cases  11  and  12,  are  associated  with 
rezoning  in  the  finite  difference  calculations.  See  Appendix 
D for  details. 

4 . 2 THEORETICAL  CONSTRAINTS  ON  THE  MULTIPOLE  COEFFICIENTS 

Our  techniques  for  synthesizing  teleseismic  body  and 
surface  wave  seismograms  require  the  multipole  coefficients 
in  the  frequency  domain.  The  appropriate  transformation  is 
given  by  (4.4).  That  is,  we  must  Fourier  transform  time 
series  like  those  in  Figures  4. 1-4. 3,  then  normalize  by  the 
Hankel  function.  Before  discussing  results  of  this  operation, 
let  us  discuss  the  constraints  on  the  form  of  the  multipole 
coefficients  that  are  imposed  by  the  theory. 

The  theoretical  constraints  on  the  solution  follow 
from  the  requirement  that  there  be  no  static  offset  in  the 
far-field.  That  is,  we  require  that 

lim  Uyp  (R, t)  = 0,  (4.6a) 

which  is  equivalent  to  requiring  that 

lim  u__  (R,w)  = ojn,  n > 0.  (4.6b) 

/\  i ** 

w-*-0 

The  far-field  portion  of  the  field  is  that  which  decays  as 
R 1.  The  relationship  between  the  far-field  displacement 
and  the  multipole  coefficients  for  the  case  being  studied 
is  given  by  (C. 47) - (C. 48)  of  Appendix  C.  We  point  out  that 
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the  conditions  (4.6)  can  be  shown  to  be  equivalent  to  re- 
quiring a finite  static  displacement  in  the  near-field. 

Comparing  (4.6b)  with  (C . 47) - (C . 4 8 ) , we  see  that  we 
must  require  that 

lim  <jjn+2,  n > 0 , (4.7) 

co-*0  m 

for  all  l and  m. 

What  are  the  implications  of  (4.7)  for  the  time  do- 
main behavior  of  the  multipole  coef f icients?  From  (4.4)  we 
observe  that 


K ,kjR» 


A ^ ^ ( R t ) ^ 


(4.8) 


where  } denotes  the  Fourier  transform  of  { }.  Since 


lim  Ap  * (R,  t)  = lim 
t-«  U) 


u l 3T t ' 


(j) 

m 

o u (— * 1 


= lim 

(j) 


(4.9) 


it  is  apparent  that  if  (R,t)  has  a finite  static  limit. 


lim^jA^  (R,t)j  - 


-1 


(4.10) 


Then,  since 


lim  — rjy 

a)-* 0 h0v  ' (k  .R) 


H+l 

0)  , 


(4.11) 


we  have 


lim  Ab(£}  (u)  = «/, 


£.m 


(4.12) 


for  those  coefficients  which  have  a finite  static  limit. 
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In  view  of  (4.7)  and  (4.12)  we  see  that  the  time  domain 
coefficients  for  l _>  2 can  have  finite  static  limits.  How- 
ever, for  the  dipole  (*-=1)  , we  must  require  that 

lim  A.(^  (R, t)  = 0,  (4.13) 

t-oo  *m 

for  the  solution  to  be  finite. 

4.3  MULTIPOLE  COEFFICIENTS  IN  THE  FREQUENCY  DOMAIN 

The  Fourier  transformation  of  terms  of  higher  order 
than  the  dipole  is  carried  out  by  assuming  that  the  value  at 
the  last  time  step  is  the  static  value.  For  the  dipole  terms, 
like  A^q  (R,t)  shown  in  Figures  4. 1-4. 3,  we  must  have  a zero 
static  value.  To  accomplish  this  the  coefficient  is  simply  set 
to  zero  for  all  times  greater  than  the  final  computed  time 
step.  For  Cases  1,  2,  3,  8,  9 and  10  this  seems  quite 
satisfactory.  The  other  six  cases  do  not  seem  to  be  approach- 
ing zero  for  long  times  and  the  approximation  is  more  severe. 
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V.  FAR-FIELD  DISPLACEMENT  SPECTRA 


For  each  of  the  fourteen  sources  listed  in  Table  3.1 
we  compute  far-field  displacement  spectra.  For  the  one- 
dimensional, spherically  symmetric  calculations  the  equivalent 
elastic  source  is  the  reduced  displacement  potential,  f (t-r/a) , 
which  is  related  to  the  displacement  by 


u (R,w) 


(5.1) 


where  all  quantities  are  Fourier  transformed.  Then  the  far- 
field  displacement  spectrum  is  defined  by 


uff(R,co) 


(5.2) 


An  entirely  analogous  procedure  is  followed  for  the 
equivalent  elastic  source  representations  for  the  cratering 
calculations.  That  is,  we  retain  only  terms  of  order  R 
in  the  expansion  as  is  explained  in  Section  C.3  of  Appendix 
C. 

The  displacement  spectra  are  presented  in  Figures 
5.1  through  5.6.  The  plots  are  log-log  in  amplitude  versus 
frequency.  Note  that  the  scale  on  the  amplitude  axis  is  in 
powers  of  10  while  the  actual  frequencies  are  printed  on  the 
abscissa.  The  spectra  are  shown  at  two  takeoff  angles  t (see 
Figure  2.1).  The  teleseismic  body  waves  are  associated  with 
takeoff  angles  near  x = 20°.  The  x = 70°  plots  are  shown  as 
being  representative  of  the  waves  trapped  as  surface  wave 
energy. 

Two  frequencies  are  singled  out  on  the  spectra  for  the 
cratering  calculations  and  are  denoted  A and  B.  The  fre- 
quency denoted  B is  associated  with  the  (approximate)  total 
time  of  the  ground  motion  data  provided  by  ATI.  That  is,  if 
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Figure  5.1.  Far-field  displacement  spectra  for  one  spherically  symmetric  and  three 
cratering  explosions  in  granite.  The  spectra  are  for  a takeoff  angle 
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Figure  5.2.  Far-field  displacement  spectra  for  one  spherically  symmetric  and  three 
cratering  explosions  in  granite.  The  spectra  are  for  a takeoff  angle 
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Figure  5.3.  Far-field  displacement  spectra  for  four  cratering  explosions  in  dry  sand- 
stone. The  spectra  are  for  a takeoff  angle  t = 20°. 
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Figure  5.4.  Far-field  displacement  spectra  for  four  cratering  explosions  in  dry 
sandstone.  The  spectra  are  for  a takeoff  angle  t = 70°. 
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Figure  5.6.  Far-field  displacement  spectra  for  one  spherically  symmetric  and  five 

cratering  explosions  in  wet  sandstone.  The  spectra  are  for  a takeoff  angle 
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the  7*u(t)  and  V*u(t)  at  a typical  station  on  the  elastic 
radius  had  their  first  non-zero  value  at  t^  = 0.4  seconds 
and  the  final  time  point  was  at  t^  = 2.5  seconds,  we  say 
that  the  frequency  B is  l/(t^-t^)  or,  for  this  example, 

0.48  Hz.  Then  B is  the  lowest  frequency  that  can  be  assoc- 
iated with  the  actual  computed  data. 

The  frequency  A is  associated  with  the  procedure  used 
to  compute  the  Fourier  transforms.  As  pointed  out  in  the  pre- 
vious section,  beyond  the  time  t^  we  assume  that  the  lowest 
order  terms  (the  dipole)  in  the  expansion  of  the  divergence 
and  curl  are  zero.  The  higher  order  terms  are  assumed  to 
remain  static  at  the  value  reached  at  the  last  time  point. 

The  amplitudes  at  all  frequencies  below  B are  dependent  on 
these  assumptions.  In  order  to  compute  the  Fourier  transform 
we  extend  the  time  histories  out  to  some  t using  the  assump- 
tions mentioned  above.  Then  A ~ l/(t  -t.). 

e l 

How  do  we  compute  values  for  frequencies  below  A?  We 
know  that  the  response  at  low  frequencies  is  dominated  by 

the  dipole  term  in  the  expansion.  Further,  in  Section  4.2 

, 2 
we  showed  that  the  dipole  term  must  behave  like  to  at  low 

frequencies  for  a bounded  solution.  This  corresponds  to  a 
flat  far-field  displacement  spectrum  at  low  frequencies.  Our 
procedure  is  then  to  extrapolate  the  amplitude  of  the  dipole 
term  from  its  value  at  A by  assuming  proportionality  to  to  . 
This  is  nearly  the  same  as  assuming  the  far-field  displace- 
ment spectrum  to  be  flat  to  long  periods  from  its  value  at 
A as  can  be  seen  in  the  figures. 

All  the  P wave  spectra  shown  are  quite  well  behaved, 
though  the  high  frequency  portion  shows  considerable  modula- 
tion. However,  some  of  the  S wave  spectra  are  rapidly 
oscillating  in  a manner  indicative  of  numerical  error.  Cases 
5,  6 and  12  are  the  worst  cases.  The  oscillations  are  ap- 
parently due  to  some  incompatibility  between  the  dipole  and 
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higher  order  terms.  We  can  ignore  this  trouble  because  only 
the  portion  of  the  S wave  spectrum  below,  say,  0.1  Hz  has  any 
influence  on  our  calculations.  The  S waves  cannot  influence 
the  body  waves  (see  Section  VI) , so  we  need  not  worry  about 
errors  in  the  higher  frequency  region. 

Finally,  in  Tables  5.1  and  5.2  we  tabulate  the  spectral 

values  that  have  the  greatest  significance  for  computing 

teleseismic  body  and  surface  waves.  The  amplitude  values  in 

these  tables  are  all  on  the  same  scale  (all  have  been  multi- 
-4 

plied  by  10  R)  and  are  best  viewed  as  relative  amplitudes. 
The  reliability  of  these  values  for  scaling  m^  and  M is 
discussed  in  Sections  VI  and  VII. 
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TABLE  5.1 

BODY  WAVE  AMPLITUDE  SPECTRA  - FAR-FIELD  P WAVE 
1 HZ  SPECTRAL  AMPLITUDES  AT  T = 20° 


Identifier 

Material 

Depth  (km) 

Spectral  Amplitude 

13 

Granite 

Spherically 

Symmetric 

1.8 

1 

II 

0.159 

1.8 

2 

N 

0.207 

2.1 

3 

It 

0.253 

1.5 

4 

Dry  Sandstone 

0.159 

0.60 

5 

It 

0.207 

0.65 

6 

It 

0.253 

0.72 

7 

Weak  Dry  Sand- 
stone 

0.207 

1.4 

14 

Wet  Sandstone 

Spherically 

Symmetric 

1.7 

8 

It 

0.053 

1.3 

9 

It 

0.159 

2.1 

10 

N 

0.207 

2.8 

11 

It 

0.253 

2.7 

12 

It 

0.531 

2.5 
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TABLE  5.2 

SURFACE  WAVE  AMPLITUDE  SPECTRA  (0.0 5 Hz) 


Identifier 

P (x=20°) 

S (t=20°) 

P (t=70°) 

S (t=70°) 

S/P  (70' 

13 

1.3 

- 

1.3 

- 

- 

1 

2.9 

2.0 

1.1 

5.4 

4.9 

2 

3.2 

5.7 

1.2 

15.6 

13.0 

3 

3.3 

0.9 

1.2 

2.7 

2.3 

4 

6.5 

1.1 

2.4 

2.8 

1.2 

5 

7.1 

5.0 

2.6 

1.4 

0.54 

6 

5.6 

5.7 

2.1 

1.6 

0.76 

7 

10.6 

0.95 

3.9 

2.4 

0.62 

14 

2.0 

- 

2.0 

- 

- 

8 

2.2 

1.7 

0.78 

4.6 

5.9 

9 

10.1 

7.4 

3.7 

20.2 

5.5 

10 

21.3 

13.6 

7.7 

37.8 

4.9 

11 

14.4 

22.0 

5.2 

60.1 

11.6 

13.7  11.8 
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VI.  BODY  WAVE  MAGNITUDE,  rr^ 

In  this  section  we  present  our  synthetic  seismograms 
for  the  fourteen  ATI  sources  and  give  the  iti^  values  for  each. 
The  synthetic  seismogram  calculations  include  the  following 
elements : 

1.  The  equivalent  elastic  sources  which  give  the 
(whole  space)  far-field  displacement  spectra 
described  in  the  previous  section  are  embedded 
in  a layered  model  of  the  crust  in  the  source 
region.  The  basic  model  used  for  the  calcula- 
tions is  tabluated  in  Table  6.1.  The  top  layer 
was  changed  to  have  the  properties  appropriate 

to  each  source  as  listed  in  Table  2.1.  In  carry- 
ing out  the  calculations  only  the  downgoing 
waves  emitted  by  the  source  are  computed;  that 
is,  no  free  surface  is  included  in  the  source 
crustal  model. 

2.  The  far-field  body  waves  emanating  from  the 
base  of  the  source  crust  and  characterized  by 
ray  parameter  p = 0.079  sec/km  are  calculated. 

3.  The  upper  mantle  is  accounted  for  by  a step 
function  response  computed  using  generalized 
ray  theory.  In  this  case  we  took  the  distance 
to  be  A = 36°  which  is  beyond  the  upper  mantle 
triplications  and  the  upper  mantle  response  is 
essentially  a constant  geometric  spreading 
factor  with  a value  of  about  0.95  x 10-4. 

4.  The  response  of  the  receiver  crustal  model  (Table 
6.2)  is  included.  In  this  case  the  receiver 
crust  has  little  effect  other  than  scaling  the 
seismogram  proportional  to  the  velocity  of  the 
top  few  kilometers. 
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TABLE  6.1 

SOURCE  REGION  CRUSTAL  MODEL 


Depth  (km)  Thickness  (km)  a (km/sec)  8 (km/sec)  p (g/cm3 ) 


* - 

1.0 

1.0 

Granite, 

Wet  or  Dry 

Sandstone 

1.7 

0.7 

4.7 

2.7 

2.6 

2.7 

1.0 

5.4 

2.8 

2.7 

. 4.0 

1.3 

5.8 

3.45 

2.8 

20.0 

16.0 

6.0 

3.50 

2.8 

TABLE  6.2 

RECEIVER 

REGION 

CRUSTAL  STRUCTURE 

Depth  (km) 

Thickness 

(km) 

a (km/sec) 

8 (km/sec) 

p (km/sec ) 

2.58 

2.58 

3.67 

2.31 

2.40 

4.84 

2.26 

5.42 

3.27 

2.60 

11.61 

6.77 

5.80 

3.45 

2.60 

20.0 

8.39 

6.00 

3.50 

2.80 
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5.  The  response  is  convolved  with  an  operator  repre- 
senting the  attenuation  and  dispersive  properties 
of  the  earth.  For  these  calculations  we  took 

t*  = 0.7  (t*  is  the  ratio  of  travel  time  to  the 
effective  path  attenuation  factor,  Q) . 

6.  The  ground  motion  is  convolved  with  the  response 
of  a standard  short  period  seismograph  system. 

The  synthetic  seismograms  are  shown  in  Figures  6.1  and 
6.2.  The  cycle  from  which  m^  is  measured  is  indicated  on 
each  record  by  a bar.  The  m^  values  and  the  period  of  the  m^ 
cycle  are  tabulated  in  Table  6.3.  The  m^  is  computed  from 

* log  | + 3.32  , (6.1) 

where  T is  the  period  tabulated,  A is  the  peak-to-peak  ampli- 
tude of  the  indicated  cycle  corrected  for  the  instrument 
response  at  the  period  T,  and  the  constant  3.32  is  the  ap- 
propriate distance  correction  factor.  Recall  that  for  all 
the  seismogram  calculations  there  is  no  free  surface  near 
the  source. 

The  seismograms  of  Figures  6. 1-6. 3 are  quite  typical 
of  short  period  teleseismic  recordings  of  explosions.  The 
only  anomaly  seems  to  be  for  the  dry  sandstone  events  for 
which  the  seismograms  are  in  Figure  6.2.  Note  that  the  ap- 
parent first  motion  in  these  seismograms  is  downward;  or  at 
least  the  upward  motion  is  too  small  to  be  noticed.  This 

seems  to  be  a consequence  of  the  negative  offset  of  the  domi- 
(4 ) 

nant  term.  A'  , that  was  discussed  in  connection  with 
1U  ' 

Figure  4.2. 

One  other  point  should  be  made  about  the  short  period 
seismograms  and  m^.  The  S wave  contribution  to  the  seismo- 
grams is  negligible.  The  only  way  the  S wave  energy  leaving 
the  source  can  affect  the  seismograms  is  through  converted 
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Figure  6.1.  Synthetic  short  period  seismograms  for  one 
spherically  symmetric  and  three  cratering 
calculations  in  granite.  The  numbers  to  the 
left  are  ground  motion  in  microns  at  1 Hz. 
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TABLE  6 . 3 


SUMMARY  OF  mb  VALUES  FOR  ATI  CRATERING  CALCULATIONS 


Identifier 

Depth 

5> 

Period 

Normalized 

(A/T)/(a*u) 

Granite 

13 

Spherically 

Symmetric 

6.29 

0.81 

0.93 

1 

0.159 

6.28 

0.78 

0.35 

2 

0.207 

6.34 

0.82 

0.82 

3 

0.253 

6.28 

0.80 

1.01 

Dry  Sandstone 

4 

0.159 

5.48 

0.96 

0.96 

5 

0.207 

5.56 

0.88 

1.09 

6 

0.253 

5.55 

0.86 

0.94 

7 (weak) 

0.207 

5.80 

0.93 

0.86 

Wet  Sandstone 

14  (1-D) 

Spherically 

Symmetric 

5.88 

0.79 

1.00 

8 

0.053 

5.66 

0.94 

0.79 

9 

0.159 

5.98 

0.92 

1.02 

10 

0.207 

6.11 

0.92 

1.05 

11 

0.253 

6.10 

1.03 

1.05 

12 

0.531 

6.15 

0.99 

1.27 
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Identifier 
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5.88 
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S'  - ' ' ^ " — | 

1 

TIME  (SEC) 

Figure  6.3.  Synthetic  short  period  seismograms  for  one 

spherically  symmetric  and  five  cratering  cal- 
culations in  wet  sandstone.  The  numbers  to 
the  left  are  ground  motion  in  microns  at  1 Hz. 
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waves  at  the  interfaces  below  the  source  (Table  6.1).  There 
is  not  even  a free  surface  to  contribute  an  sP  phase.  For 
the  velocity  contrasts  in  the  structure  and  the  steep  take- 
off angles  for  teleseismic  body  waves,  very  little  of  the  S 
wave  energy  is  partitioned  into  early  arriving  P waves. 

It  is  interesting  to  examine  how  the  scales  with 
the  far-field  displacement  spectra  discussed  in  Section  V. 

To  first  order  we  expect  to  have  [Bache,  et  al . , 1976], 

= log  uj , (6.2) 

where  u is  the  displacement  spectrum  at  the  controlling  fre- 
quency and  a is  the  P wave  velocity  at  the  source.  Using 
s ^ 

the  values  of  from  Table  2.1  and  u (1  Hz)  from  Table  5.1, 
we  normalize  the  amplitudes  and  scale  to  the  value  for 
calculation  14.  The  normalized  scaled  amplitudes  are  tabu- 
lated in  the  last  column  of  the  table.  We  see  that  all 
values  are  within  + 27  percent  of  the  normalized  amplitude 
for  calculation  14.  Further,  eleven  of  the  fourteen  are 
scaled  to  within  15  percent  by  (6.2).  The  discrepancies  are 
mainly  attributed  to  the  fact  that  we  are  normalizing  to  the 
spectral  amplitude  at  1 Hz  while  the  frequencies  controlling 
the  seismogram  amplitudes  range  from  0.97  - 1.32  Hz.  From 
Figures  4.1  - 4.6  we  see  that  the  P wave  spectra  are  complex 
and  rapidly  changing  in  this  region. 
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VII.  SURFACE  WAVE  MAGNITUDE,  M 

s 

In  this  section  we  present  our  synthetic  long  period 
seismograms  for  the  fourteen  ATI  sources  and  give  the  Mg 
values  for  each.  The  synthetic  seismogram  calculations  in- 
clude the  following  elements: 

1.  The  method  for  the  surface  wave  calculations  is 
described  by  Harkrider  [1964]  and  Harkrider  and 
Archambeau  [1977] . The  same  equivalent  elastic 
source  formulation  used  for  the  body  waves  is 
used  for  the  surface  waves.  Once  again,  only 
the  down-going  waves  from  the  source  are  in- 
cluded in  the  calculations. 

2.  Two  crustal  models  are  used  for  the  path,  one 
for  the  very  near  source  region  and  one  for  the 
remainder  of  the  path  to  the  receiver.  The 
average  path  model  is  one  proposed  for  North 
America  by  McEvilly  [1964].  The  only  difference 
between  the  models  is  that  the  top  three  kilo- 
meters of  the  source  region  crustal  model  is  re- 
placed by  the  ATI  granite,  wet  or  dry  sandstone, 
depending  on  the  source  material.  For  the  long 
periods  controlling  teleseismic  Mg  the  reflection 
coefficient  for  Rayleigh  waves  passing  across  the 
boundary  between  the  source  and  average  path 
crustal  model  is  close  to  unity. 

3.  The  ground  motion  is  convolved  with  the  response 
of  an  LRSM  long  period  seismometer.  A Q operator 
which  has  only  a minor  effect  is  also  included. 
The  Q model  is  that  of  Tryggvason  [1965] . 

4.  The  seismograms  were  synthesized  at  a range  of 
3000  km.  The  Mg  was  computed  using  the  formulas 
of  Marshall  and  Basham  [1972] . For  this  range 
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the  formula  reduces  to 

M = log  A + 1. 38  + P (T) 
s 

where  A is  the  maximum  amplitude  (zero-to-peak) 
of  the  signal  with  period  near  20  seconds  and 
P (T)  is  a period  dependent  path  correction  tabu- 
lated by  Marshall  and  Basham  [1972] . The  cor- 
rection is  quite  small  for  periods  near  20 
seconds . 

The  vertical  component  Rayleigh  wave  synthetic  seismo- 
grams are  shown  in  Figures  7.1  - 7.3.  Two  seismograms  are 
shown  for  the  spherically  symmetric  contained  explosions  in 
granite  and  wet  sandstone  (Cases  13  and  14).  In  one  of  these 
(13b,  14b)  only  the  downward  waves  are  included  in  the  cal- 
culation. The  seismograms  are  then  the  analog  of  the  Cases 
13  and  14  body  wave  seismograms  of  Figures  6.1  and  6.3. 

Since  the  Mg  is  relatively  insensitive  to  depth  for  contained 
explosions  (unlike  m^  which  is  strongly  affected  by  the  pP 
phase) , we  also  compute  seismograms  for  fully  contained  explo- 
sions at  a depth  of  200  meters.  The  seismograms  are  given 
as  Cases  13a  and  14a. 

The  important  data  from  the  seismograms  of  Figures 
7.1  - 7.3  are  summarized  in  Table  7.1.  The  phase  at  which 
the  amplitude  for  Ms  was  measured  is  indicated  on  the  seis- 
mograms by  a bar.  The  period  of  this  phase  is  given  in  the 
table.  Also  listed  in  the  table  is  the  spectral  amplitude 
of  the  true  ground  motion  at  a period  of  20  seconds.  The 

A 

difference  between  Mg  and  log  A indicates  the  consistency  to 

which  the  Mg  measurement  represents  a true  measurement  of  the 

energy  at  frequencies  in  this  range.  These  differences  are 

tabulated  in  the  last  column  of  Table  7.1.  We  see  _hat  the 

Mg  values  are  quite  consistent  with  the  spectral  measurements 

with  the  spread  between  the  maximum  and  minimum  values  being 

0.14  M units, 
s 
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Figure  7.1.  Vertical  component  Rayleigh  wave  seismograms 

for  a spherically  symmetric  and  three  cratering 
explosions  in  granite.  The  numbers  at  the  left 
are  displacements  in  microns  at  25  seconds. 
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Identifier 


M 


s 


TIME  (SEC) 


Figure  7.2.  Vertical  component  Rayleigh  wave  seismograms 

for  four  cratering  explosions  in  dry  sandstone. 
The  numbers  at  the  left  are  displacements  in 
microns  at  25  seconds. 
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Vertical  component  Rayleigh  wave  seismograms 
for  a spherically  symmetric  and  five  cratering 
explosions  in  wet  sandstone.  The  numbers  at  the 
left  are  displacements  in  microns  at  25  seconds. 
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TABLE  7.1 

SUMMARY  OF  M VALUES  FOR  ATI  CRATERING  CALCULATIONS 
s 


Identifier 

Depth  (km) 

M 

s 

M Period 
s 

Log  Spectral 
Amplitude, 
log  A (20  sec) 

Mg-log 

Granite 

13A 

1-D  Total  Source 

4.57 

20.1 

1.92 

2.65 

13B 

1-D  Down  Waves 
Only 

4.34 

21.5 

1.62 

2.72 

1 

0.159 

4.01 

20.0 

1.30 

2.71 

2 

0.207 

4.95 

23.6 

2.25 

2.70 

3 

0.253 

4.72 

21.3 

2.03 

2.69 

Dry  Sandstone 

4 

0.159 

4.66 

22.0 

1.92 

2.74 

5 

0.207 

4.68 

22.7 

1.92 

2.76 

6 

0.253 

4.60 

22.5 

1.83 

2.77 

7 ('weak) 

0.207 

4.83 

22.0 

2.07 

2.76 

Wet  Sandstone 

14A 

1-D  Total  Source 

4.17 

22.8 

1.52 

2.65 

14B 

1-D  Down  Waves 
Only 

4.26 

20.0 

1.53 

2.73 

8 

0.053 

4.19 

21.2 

1.43 

2.76 

9 

0.159 

4.71 

21.2 

2.02 

2.69 

10 

0.207 

4.99 

19.3 

2.35 

2.64 

11 

0.253 

5.07 

20.1 

2.44 

2.63 

12 

0.531 

5.00 

19.8 

2.34 

2.66 
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Which  portions  of  the  wave  field  are  most  important 
for  the  Ms?  This  is  an  important  question  since  we  have  more 
confidence  in  some  portions  of  the  calculation  than  others. 

In  particular,  we  expect  the  P waves  to  be  computed  much 
more  accurately  than  the  S waves  in  our  equivalent  elastic 
source  calculation. 

In  Figure  7.4  we  show  surface  wave  seismograms  for 
several  cases  in  which  the  S wave  portion  of  the  field  has 
been  suppressed.  The  important  data  from  these  seismograms 
is  tabulated  in  Table  7.2.  We  see  that  the  relationship 
between  spectral  amplitudes  and  Ms  is  essentially  the  same 
as  for  the  total  field  seismograms  in  Table  7.1. 

The  next  step  is  to  compare  the  S-suppressed  seismo- 
grams to  the  total  wave  records  and  relate  the  results  to 
the  amount  of  S waves  present  in  the  source.  This  is  done 
in  Table  7.3.  The  S/P  ratios  are  for  the  spectra  at  t = 70° 
and  are  taken  from  Table  5.2.  The  cases  we  have  selected 
span  the  range  of  possibilities.  For  granite  and  wet  sand- 
stone we  have  the  cases  with  the  most  and  least  S wave 
component.  The  dry  sandstone  case  has  the  least  S waves  of 
all  those  studied. 

From  Table  7.3  we  see  that  if  the  S waves  are  over- 
estimated for  Case  2,  and  study  of  the  spectra  in  Figure  4.2 
indicates  they  may  well  be,  we  may  be  overestimating  M by 
0.2-0. 3 magnitude  units.  Case  2 is  by  far  the  worst  case. 

The  S waves  seem  to  have  little  effect  on  M for  the  dry 

s 

sandstone  calculations.  For  wet  sandstone  the  smaller  S/P 
ratio  acts  to  degrade  Mg  while  the  larger  ratio  acts  to 
increase  it.  For  the  other  three  cases  with  S/P  between  the 
extremes,  the  S waves  are  not  likely  to  be  too  important. 
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TABLE  7 . 2 

SURFACE  WAVE  DATA  FROM  THE  SEISMOGRAMS  OF  FIGURE  7.4 


Identifier 


P . , 
M Period 
s 


Log  Spectral 
Amplitude, 
log  A (20  sec) 


M^-log  Ap 


Granite 


2 

4.55 

21.6 

1.89 

2.66 

3 

4.56 

21.7 

1.96 

2.66 

Dry  Sandstone 

5 

4.64 

21.8 

1.93 

2.71 

Wet  Sandstone 

10 

5.12 

22.1 

2.41 

2.71 

11 

4.96 

22.3 

2.24 

2.72 

TABLE 

7.3 

COMPARISON  OF  DATA  FOR  TOTAL  SOURCE  AND  S- 
SURFACE  WAVE  SEISMOGRAMS 

-SUPPRESSED 

Identifier 

S/P 

T P 

M-M  . 

s s log 

■'T  ~P 

A -log  A 

Granite 

2 

13.0 

0.40 

0.36 

3 

2.3 

0.16 

0.13 

Dry  Sandstone 

5 

0.54 

0.04 

1 

o 

• 

o 

H* 

Wet  Sandstone 

10 

4.9 

-0.13 

-0.06 

11 

11.6 

0.11 

0.20 
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Identifier 

M 

s 

i 

2 

4.95 

•4 

Total  Source 

' 

7 

2 

- - A 

A A 'A 

1,. 

4.55 

1 

N 

a 

P Waves  Only 

I 

V 

i 

i 

3 

_ A 

ft  % 

4.72 

Total  Source 

* 

y 

j 

, 

3 

_ ^ A 

A \ l\ 

4.56 

>■* 

P Waves  Only 

J Jij 
1 

j " 
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Figure  7.4.  Vertical  component  Rayleigh  wave  seismograms  at 
3000  km  for  five  of  the  cratering  calculations. 
The  first  record  of  each  pair  is  taken  from 
Figures  7. 1-7. 3.  The  second  record  is  for  the 
same  source  with  the  S wave  portion  suppressed. 
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VIII.  CONTRIBUTION  OF  EJECTA  FALLBACK  TO  TELESEISMIC 

MAGNITUDES 


8 . 1 BOUNDARY  CONDITIONS 

In  carrying  out  the  cratering  calculations,  ATI  also 
kept  track  of  the  material  ejected  from  the  crater.  This 
material  falls  back  to  earth  and  the  impact  is  a source  of 
seismic  waves.  We  computed  the  teleseismic  body  and  surface 
waves  for  the  ejecta  fallback  contribution  to  the  cratering 
source  and  the  results  are  discussed  in  this  section. 

For  the  ejecta  field  ATI  provided  stress  histories 
on  the  free  surface.  At  each  point  (or  ring  for  this  axisym- 
metric  problem)  on  the  free  surface,  r = r' , the  time  history 
of  the  normal  stress  could  be  expressed  as 


N 


-E 

i=l 


A. 

i 


(Ti) 


where 


t . 

l 


9 


( i#  M < \ • 

n (x)  = ) 

( o,  I x I > J , 


(8.1) 


and  At  is  the  time  step,  while  A^  is  the  amplitude  of  the 
normal  stress  at  the  ifch  time  step. 

We  observe  that 


JF{n(Tt>>  -/  n(llii) 


,ooAt, 

i I .-1*  dt  . 2 -n  ‘T-’  ~luti 

CO 


(8.2) 
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Then  the  normal  stress  boundary  condition  applied  at  r = r'  is 


PQ  (r  ' , co)  = 


2 sin  (2§£)  N 


i oj  t . 

- . 


w 


(8.3a) 


i=l 


and  for  the  shear  stress  it  is 


Q0  (r',co) 


0 „ , coAt.  N 

2 sin  ) 


co 


-Ebi 


— iuj t . 
1 


(8.3b) 


i=l 


where  B.  is  the  amplitude  of  the  applied  shear  stress  at  the 
ifc^  time  step. 

The  boundary  conditions  (8.3)  are  in  the  correct 
form  for  applying  the  theory  developed  in  Appendix  E.  The 
total  stresses  on  the  surface  enter  the  solution  via  the 
spatial  integrals  (33)  of  that  Appendix.  The  quadrature  is 
straightforward  since  the  spatial  dependence  of  the  stresses 
is  also  in  terms  of  rectangle  functions.  For  example,  for 
the  normal  stresses  it  is  necessary  to  compute 


<“k> 


w 

-f 


po(r 


,ui) 


r ')  r'dr' 


(8.4) 


for  each  frequency  for  subsequent  transformation  back  to 
the  time  domain.  Here  kR^  represents  the  wave  number  for 
the  fundamental  mode  at  frequency  Then 


p(“k> 


b N 


-f  Yk  Vr'>  e 


-iu),  t . 
k l 


a i=l 


V\  r‘> 


N . . b 

lw,  t . 

Yk  7 e 


i=l  a K 


r 'dr  ' , 

(8.5) 
r 'dr  ' , 
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where 


Yk  = 


At  GJ 

2 sin  ( — 

GL 


Since  A^(r')  is  piecewise  constant: 


Ai(r')  = A^  j , a ^ < r ' <_  b^  , 


this  reduces  to 


N . M . 


i=l 


Pi  A 


b.\  " a.J  fkn  a . \ 

3 M \ V 3 1\  \ 3) 


(8.6) 


A similar  procedure  is  followed  for  the  Q(gj,  ) except  that 

K 

quadratures  are  used  to  evaluate  the  integrals 
b . 


/ 


(kR  r')  r'dr' . 

aj  j 


Then  with  P(gj^)  and  Qtw^)  we  can  proceed  to  compute  tele- 
seismic  body  or  surface  waves  using  the  formulas  (35)  and 
(44)  of  Appendix  E. 


8 . 2 SURFACE  WAVES  FROM  EJECTA 

Surface  wave  seismograms  were  computed  for  the  ejecta 
field  for  several  cases.  The  seismograms  are  shown  in  Figure 
8.1.  The  seismograms  were  computed  at  a range  of  4000  km 
rather  than  the  3000  km  range  used  for  the  records  of 
Section  VII.  In  order  to  see  the  effect  of  range  on  M 
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Figure  8.1.  Vertical  component  Rayleigh  wave  seismograms  due 
to  the  ejecta  fallback  from  six  cratering  explo- 
sions. The  numbers  at  the  left  are  displacements 
in  microns  at  25  seconds. 
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(the  distance  correction  factor  in  the  Mg  formula  may  differ 
from  the  actual  amplitude  decay) , we  computed  a seismogram 
at  3000  km  for  two  of  the  cases. 

The  data  for  the  ejecta  surface  wave  seismograms  are 
summarized  in  Table  8.1.  The  explosion  Mg  is  from  Table 
7.2.  The  difference  between  the  two  M values  gives  a rea- 
sonable  approximation  to  the  relative  amplitudes.  For 
example,  when  the  difference  is  2,  the  amplitude  of  the  ejecta 
surface  wave  at  periods  near  20  seconds  is  about  two  orders 
of  magnitude  smaller  than  that  of  the  surface  wave  from  the 
explosion  coupling  directly  into  the  ground.  In  fact,  for 
the  two  cases  (1  and  3)  for  which  we  have  seismograms  at  the 
same  range  (3000  km)  , we  compare  the  amplitudes  directly  in 
Table  8.2. 

It  will  suffice  to  examine  Case  1 since  for  this  case 
we  have  the  largest  ejecta  surface  wave  together  with  the 
smallest  direct  explosion  surface  wave.  The  comparison  of 
M values  in  Table  8.1  suggests  that  the  amplitudes  of  the 
20  second  portion  of  the  wavetrain  should  differ  by  a fac- 
tor of  about  6.5.  The  actual  ratio  is  available  from  the 
numbers  in  Table  8.2  and  is  about  6.8.  The  maximum  ampli- 
tudes differ  by  considerably  more,  however,  with  the  ratio 
being  about  34. 

From  the  analysis  of  the  previous  paragraph  we  see 
that  if  the  ejecta  contribution  were  exactly  in  phase  with 
the  direct  explosion  surface  wave,  an  unlikely  circumstance, 
the  maximum  contribution  (Case  1)  would  change  the  ampli- 
tude by  15  percent  or  0.06  magnitude  units.  For  all 
the  other  cases  the  possible  ejecta  contribution  is  at 
least  an  order  of  magnitude  less  than  this.  Therefore,  we 
feel  justified  in  ignoring  the  ejecta  contribution  to  the 
surface  waves. 

I 
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TABLE  8.1 

SURFACE  WAVE  DATA  FROM  EJECTA  SEISMOGRAMS 


Identifier 

Ejecta  Mg 
Range  (km) 

Ejecta 

Ms 

Explosion 

Ms  Difference 

Granite 

1 

4000 

3.22 

4.01 

0.79 

1 

3000 

3.20 

4.01 

0.81 

2 

3000 

3.21 

4.95 

1.74 

3 

4000 

2.90 

4.72 

1.82 

3 

3000 

2.66 

4.72 

2.06 

Dry  Sandstone 

4 

4000 

2.17 

4.66 

2.49 

5 

4000 

1.51 

4.68 

3.17 

Wet  Sandstone 

9 

4000 

2.91 

4.71 

1.80 

TABLE  8. 

2 

COMPARISON  OF  SURFACE  WAVE  AMPLITUDES  AT  3000  km 

(Amplitudes  are  in  nm  and 

are  corrected  for  instrument  res- 

ponse ) 

M 

s 

M 

s 

Maximum 

Period  of 

Identifier 

Amplitude 

Period 

Amplitude 

Maximum  Cycle 

1 Explosion 

425.0 

20.0 

5435.0 

11.3 

1 Ejecta 

62.0 

22.0 

158.0 

16.3 

3 Explosion 

2081.0 

21.8 

11,242.0 

14.3 

3 Ejecta 

17.7 

22.0 

70.0 

15.6 
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8 . 3 BODY  WAVES  FROM  EJECTA 

We  saw  in  the  previous  paragraph  that  the  ejecta  is 
unlikely  to  have  any  noticeable  effect  on  the  surface  waves 

A A 

or  Mg.  Since  the  source  terms,  P(co)  and  Q(w)  of  Appendix  E, 
are  the  same  for  both  cases,  they  would  have  to  be  strongly 
peaked  (values  at  1 Hz  larger  than  at  low  frequencies)  for 
there  to  be  any  chance  of  an  effect  on  the  short  period  body 
wave  records.  However,  as  is  reasonably  clear  from  (8.3)  or 
(8.6),  the  source  spectra  have  their  largest  values  at  low 
frequencies  and  decay  rapidly  with  increasing  frequency.  The 
spectral  values  at  the  m^  frequencies  are  down  by  at  least 
an  order  of  magnitude  from  the  values  at  Mg  frequencies. 
Therefore,  we  can  ignore  the  ejecta  contribution  to  the  short 
period  recordings. 
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IX.  TELESEISMIC  MAGNITUDES  FOR  37.5  AND  600  kt 
9.1  SOURCE  SCALING 

If  we  assume  cube-root  scaling  is  applicable,  we  can 
scale  our  results  to  obtain  and  Mg  for  cratering  explo- 
sions with  yields  different  than  the  150  kt  at  which  the 
original  calculations  were  made.  The  amplitudes  of  the  curl 
and  divergence  of  the  displacement  field  are  unchanged  by 
the  scaling,  but  the  distance  and  time  scales  are  multipled 
by  a factor  (Wn/150)^^,  where  is  the  new  yield. 

The  basic  representation  of  the  source  for  computing 
teleseismic  body  and  surface  waves  is  in  terms  of  the  multi- 
pole coefficients,  (cj)  , in  the  frequency  domain.  It  is 

easily  shown  that  cube-root  scaling  implies  that  the  scaled 
multipole  coefficients  for  the  new  yield  are  related  to  the 
old  coefficients  by 


where 


To  compute  seismograms  for  the  scaled  source  we  need  to 
scale  the  multipole  coefficients  as  indicated  and  then  to 
follow  through  the  procedures  outlined  in  Sections  VI  and 
VII. 


9.2  BODY  WAVES  FOR  SCALED  SOURCES 

We  compute  seismograms  for  the  sources  of  Table  3.1 
at  two  new  yields,  37.5  and  600  kt.  The  scale  factor  is 
then  the  cube-root  of  4 and  its  reciprocal.  The  seismograms 
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Figure  9.1. 


Synthetic  short  period  seismograms  for  the 
granite  sources  scaled  to  37.5  kt. 
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Figure  9.2. 


Synthetic  short  period  seismograms  for  the 
granite  sources  scaled  to  600  kt. 
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Figure  9.3. 


Synthetic  short  period  seismograms  for  the  dry 
sandstone  sources  scaled  to  37.5  kt. 
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Figure  9.4.  Synthetic  short  period  seismograms  for  the  dry 
sandstone  sources  scaled  to  600  kt. 
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Figure  9.5. 


Synthetic  short  period  seismograms  for  the  wet 
sandstone  sources  scaled  to  37.5  kt. 
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Figure  9.6.  Synthetic  short  period  seismograms  for  the  wet 
sandstone  sources  scaled  to  600  kt. 
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TABLE  9.1 


SUMMARY  OF  DATA  FOR  37.5  AND  600  kt  EXPLOSIONS 


37.5  kt 

600  kt 

Identifier 

% 

Period 

Period 

Granite 

13  (1-D) 

5.70 

0.79 

6.82 

0.85 

1 

5.54 

0.67 

6.82 

0.91 

2 

5.68 

0.73 

6.82 

0.92 

3 

5.54 

0.69 

6.80 

0.92 

Dry  Sandstone 

4 

4.94 

0.79 

5.65 

0.79 

5 

5.04 

0.78 

5.70 

0.78 

6 

5.07 

0.79 

5.69 

0.79 

7 (weak) 

5.31 

0.80 

6.13 

0.77 

Wet  Sandstone 

14  (1-D) 

5.27 

0.79 

6.44 

0.80 

8 

5.00 

0.72 

6.06 

0.77 

9 

5.43 

0.79 

6.36 

0.77 

10 

5.68 

0.83 

6.50 

0.78 

11 

5.69 

0.84 

6.18 

0.71 

12 

5.62 

0.79 

6.27 

0.87 
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are  shown  in  Figures  9. 1-9. 6.  The  data  from  these  seismo- 
grams are  summarized  in  Table  9.1.  Once  again,  the  phase  at 
which  the  m^  is  measured  is  indicated  on  the  seismograms. 

In  Table  9.2  we  compare  the  data  for  the  three  yields. 
We  take  the  difference  between  the  m^  for  150  kt  and  the  m^ 
values  for  37.5  and  600  kt.  The  column  marked  "Factor"  is 
then  the  antilog  of  this  difference.  We  also  show  the  ratios 
of  the  periods  of  the  m^  phases  for  the  two  yield  pairings. 

If  the  source  spectra  were  flat  throughout  the 
band  of  interest,  the  "Factor"  values  in  Table  9.2  would 
all  be  4.0.  The  values  that  appear  in  the  table  indicate 
the  complexity  of  the  source.  This  is  clearly  illustrated 
by  comparing  the  seismograms  for  the  same  cases  in  Figures 
6. 1-6. 3 and  Figures  9. 1-9. 6. 

Comparing  the  37.5  kt  records  to  the  150  kt  records, 
the  behavior  is  more-or-less  as  expected.  Sometimes  the 
amplitude  scales  by  more  than  a factor  of  4,  sometimes  by 
less.  The  periods  are  shorter  on  the  37.5  kt  records  by  a 
fairly  uniform  ratio  throughout.  For  the  600  kt  granite 
events,  the  results  are  also  within  the  range  of  expected 
behavior . 

The  comparison  of  600  kt  and  150  kt  records  for  the 

wet  and  dry  sandstone  gives  some  surprising  results.  The 

amplitude  ratios  are  considerably  less  than  four  for  all 

the  cratering  events.  Further,  the  periods  for  the  600  kt 

m^  are  actually  shorter  than  for  the  150  kt  m^.  In  fact, 

the  periods  for  the  37.5  kt  and  600  kt  sandstone  seismograms 

are  about  the  same.  The  anomaly  seems  to  be  due  to  the 

interference  pattern  in  the  mfe  portion  of  the  wavetrain. 

A strong  interfering  phase  is  clearly  visable  on  the  600  kt 

records  for  Cases  4,  5,  6,  11  and  12.  Looking  at  the  data 

of  Table  9.2,  these  are  the  events  for  which  the  m,  is 

b 

smallest  compared  to  the  expected  values  from  scaling  the 
150  kt  amplitudes  by  four. 
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The  failure  of  the  body  wave  amplitudes  to  scale 
exactly  with  yield  is  a function  of  the  complexity  of  the 
source  spectrum  in  the  short  period  region.  However,  the 
long  period  or  Mg  portion  of  the  spectrum  is  flat  (Figures 
5. 1-5. 6).  Therefore,  there  is  no  reason  to  expect  the  Ms 
values  to  change  by  any  value  other  than  the  log  of  the 
yield  ratio  or  + 0.60. 
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APPENDIX  A 

EQUIVALENT  ELASTIC  SOURCE 


The  radiation  field  exterior  to  any  kind  of  volume 
source  in  a homogeneous  medium  can  be  represented  in  terms  of 
an  expansion  in  spherical  harmonics.  Archambeau  (1968)  seems 
to  have  been  the  first  to  recognize  the  usefulness  of  this 
fundamental  result  and  to  apply  it  to  geophysical  problems. 
The  expansion  in  spherical  harmonics  gives  a compact  equiva- 
lent elastic  source  representation  of  quite  general  character 
and  nearly  any  proposed  seismic  source  model  can  be  cast  in 
this  form.  A brief  description  of  this  source  representation 
and  its  compatibility  with  commonly  used  source  theories  is 
the  subject  of  this  section. 

The  Fourier  transformed  equations  of  motion  in  a homo- 
geneous, isotropic,  linearly  elastic  medium  may  be  written 


where  u is  particle  displacement  and  k^  and  kg  are  the 

conpressional  and  shear  wave  numbers.  The  Cartesian  poten- 
— (4 ) — 

tials  x and  x are  defined  by 


X u , 


(2) 


and  may  be  easily  shown  to  satisfy  the  wave  equation 
V2X(j)  + k2X(j)  = 0,  j = 1,  2,  3,  4, 


(3) 
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where  k = k = co/a  and  k.  = k0  = oo/B  for  i = 1,2,3.  This 

4 Ufc  1 p 

equation  has  as  a solution  the  following  expansion  in  spherical 
eigenfunctions  (e.g.,  Morse  and  Feshbach  (1953)), 


X(j)  <R,w)  = ^2  h{z2)  (k j R) 

2=0 


cos  mi? 


+ B^  (w)  sin  m<J)  P™(cos6)  , (4) 

(2) 

where  the  h,  are  spherical  Hankel  functions  of  the  second 
m 

kind  and  the  P , are  associated  Legendre  functions.  The 

vector  R has  as  components  the  spherical  coordinates  R,9,<j>. 

Equations  (4),  together  with  (1),  provide  an  elastic 

point  source  representation  of  the  (outgoing)  displacement  field. 

The  values  of  the  multipole  coefficients,  A^j^  (u>)  , (cj)  , 

j = 1,2, 3, 4,  prescribe  the  displacement  field  at  all  points  in 

the  homogeneous  medium  where  (1)  applies.  This  point  source 

representation  can  be  viewed  as  a generalized  form  for  a sum 

of  a monopole  or  center  of  dilatation  (2  = 0),  a dipole  or 

couple  (2  = 1)  , a quadrupole  or  double-couple  (2.  = 2),  etc. 

For  example,  a center  of  dilatation  is  represented  by  a single 
( 4 ) 

coefficient  A , while  for  a horizontal  double-couple  the 
o o 

nonzero  coefficients  are  -A  ^ = a^  and  B^. 

21  21  22  22 

Description  of  the  character  of  the  elastic  field  gene- 
rated by  seismic  sources  is,  of  course,  a basic  geophysical 
problem.  For  this  paper  it  is  convenient  to  discuss  seismic 
source  descriptions  in  three  categories: 
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1.  Those  obtained  using  finite  difference/finite 
element  numerical  methods. 

2.  Analytical  source  models  of  relaxation  type. 

3.  Dislocation  source  models. 

With  numerical  methods  one  can  attempt  to  directly  in- 
clude complexities  of  the  source  mechanism  in  a deterministic 
computational  scheme.  For  example,  finite  difference  methods 
have  been  extensively  used  to  compute  the  propagating  shock 
wave  due  to  an  underground  nuclear  explosion  (e.g..  Cherry, 
et  al . , 1974)  . In  this  case  the  nonlinear  behavior  of  the 
rock  under  high  stress  loading  determines  the  character  of 
the  seismic  signal.  If  the  source  region  can  be  assumed  to 
be  embedded  in  a medium  in  which  (1)  applies,  an  equivalent 
elastic  source  of  the  form  (4)  can  be  obtained  from  the  out- 
going displacement  field.  This  is  indicated  schematically  in 
Figure  1.  Briefly,  the  procedure  is  to  monitor  the  outgoing 
displacement  field  or,  alternatively,  the  potentials,  x ^ , 

A 

on  a spherical  surface  of  radius  R.  Using  the  orthogonality 
of  the  spherical  harmonics,  these  potentials  are  related  to 
the  multipole  coefficients  by 


<->) 

I Bim  <“> 

where 


')lm 


2 TT  7T 


'im 


,m 


'C2 ~ J J X(^(R,w)  P‘“(cos9)  { } sin9d6d<(> , 


h£  (kjR)  •'o  "0 


j cos  mg ) 

/ r*  «»  ■ 

| sin  mgj 


(5) 


(2&+1) (&-m) 1 
2 tt  ( £+m) 

= (2  £+1 ) /4r  . 


, m ^ 0 , 
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Use  of  this  procedure  for  linking  nonlinear  finite  dif- 
ference source  calculations  with  analytical  wave  propagation 
techniques  was  suggested  to  the  first  author  by  Archambeau 
(1973,  personal  communication),  and  has  since  been  implemented 
for  a number  of  complex  explosion  geometries  (Cherry,  et  al . , 
1975,  1976a),  and  for  a three-dimensional  finite  difference 
simulation  of  stick-slip  earthquake  faulting  (Cherry, 
et  al . , 1976b) . The  number  of  terms  required  for  the  expan- 
sion (4)  to  converge  depends  on  the  symmetry  of  the  source 
radiation  at  frequencies  of  interest.  The  most  elementary  ap- 
plication of  the  method  is  for  one-dimensional  (spherically 
symmetric)  explosion  source  calculations.  For  such  problems 
the  elastic  field  is  often  described  by  a reduced  displacement 
potential  defined  by 


u(R,t>  = 


V ( t-R/a) 


(6) 


Applying  the  Fourier  transform  and  comparing  to  (1)  together 
with  (4),  it  is  easily  derived  that 


a(4)m  = 
0 0 


-i  k!  Y(u)  , 
a 


(7] 


which  shows  the  equivalence  between  the  reduced  displacement 
potential  and  the  monopole.  For  more  complex  sources  such 
as  an  explosion  in  an  axisymmetric  tunnel  (Cherry,  et^  al , , 
1975)  or  several  explosions  detonated  simultaneously  (Cherry, 
et  al . , 1976a)  quadrupole  and  higher  order  terms  occur  in 
the  expansion.  When  an  earthquake  source  is  computed,  the 
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leading  term  is,  as  expected,  the  quadrupole  (Cherry,  et  al., 


1976b) . 
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APPENDIX  B 

APPLICATION  OF  THE  MULTIPOLAR  EXPANSION  TECHNIQUE 
TO  CRATERING  EXPLOSIONS 


In  Appendix  A is  outlined  the  procedure  for  obtaining 
an  equivalent  elastic  source  representation  of  an  arbitrary 
volume  source.  In  particular,  we  evaluate  the  integrals 

(A. 5)  to  obtain  the  multipolar  coefficients  (A,f^  (w)  , 

( 1 ) ?'m 

B„m  (w) ) that  define  the  equivalent  elastic  source.  With 

these  multipoles  the  displacement  can  be  computed  from 

( A . 1 ) and  (A. 4)  at  any  point  in  space. 

The  procedure  outlined  in  Appendix  A gives  a unique 
and  exact  representation  of  the  outgoing  wave  field  when 
the  following  conditions  are  satisfied: 

1.  The  spherical  surface  at  the  "elastic  radius" 
is  in  a homogeneous,  isotropic,  linearly 
elastic  medium. 

2.  No  energy  travels  inward  through  this  surface. 

3.  A sufficient  number  of  terms  are  computed  to 
insure  convergence  of  the  infinite  series,  Eq. 

(A. 4)  of  Appendix  A. 

The  ideal  conditions  listed  above  are  certainly  not 
satisfied  in  the  problem  of  interest.  The  most  obvious 
and  important  difficulty  is  the  presence  of  a free  surface 
To  understand  the  effect  of  the  free  surface,  we  should 
first  describe-  the  relevant  features  of  how  our  technique 
works.  First,  the  multipolar  form  can  only  represent  waves 
that  are  solutions  to  the  whole  space  wave  equation;  that 
is,  body  waves  of  compressional  or  shear  type.  Due  to  the 
presence  of  the  free  surface,  there  will  be  surface  waves 
of  the  Rayleigh  type  in  the  numerical  solution.  These  waves 
can  never  be  properly  represented  in  our  multipolar  solu- 
tion. The  solution  will  "see"  the  Rayleigh  wave  displacements 
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as  some  combination  of  P and  S waves  and  will  try  to  repre- 
sent them  in  that  form.  In  subsequent  paragraphs  we  will 
explore  the  effect  of  the  Rayleigh  waves  in  more  detail  and 
argue  that  they  do  not  seriously  compromise  our  solution. 


The  second  characteristic  of  our  equivalent  elastic 
source  determination  technique  that  needs  to  be  understood 
is  the  way  in  which  energy  is  partitioned  into  the  various 
multipoles.  It  is  the  symmetry  of  the  outgoing  displacement 
field  that  controls  this  partitioning  via  the  integral 
(A. 5).  We  see  from  Eq.  (A. 5)  that  we  operate  on  the  diver- 
gence (x^(RfW))  and  curl  (x^(R,w),  i = 1,2,3)  of  the 
displacement  field  separately.  As  an  illustrative  example, 
let  us  assume  that  we  have  an  axisymmetric  source  in  a ho- 
mogeneous whole  space  and  examine  how  the  propagating  di- 
vergence is  partitioned  into  the  multipole  coefficients 


— (4 ) 

Since  the  x is  independent  of  <p , the 
(4) 


only  non-zero  multipoles  are  A„  ' (w) . Then,  from  (A. 4), 
the  divergence  is  written  as  follows: 


X(4)(R,w)  = ^ hj£2)  (kaR)  AZ0)(UJ)  p°(cos9>  (B.l) 

£=0 

That  is,  each  frequency  component  of  the  source  generated 

displacement  field  is  made  up  of  a sum  of  terms  with  symmetry 

specified  by  the  Legendre  polynomials.  The  radiation  pattern 

for  the  first  several  terms  is  shown  in  Figure  B.l.  Then, 

(4 ) 

for  example,  if  the  divergence,  X , of  the  displacement 
field  exhibits  considerable  spherical  symmetry,  the  solution 
will  be  dominated  by  the  first  term  in  B.l  and  the  higher 
order  terms  will  represent  the  deviation  from  this  symmetry. 
The  same  kind  of  process  operates  for  the  curl  and  the  argu- 
ments are  easily  extended  to  a more  general  (not  axisymmetric) 
source . 
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Figure  B.l.  The  radiation  patterns  for  the  first  several 
terms  of  the  expansion  of  the  divergence  for 
an  axisymmetric  source.  The  p£  is  the  radia- 
tion pattern  for  the  aI4) (w)  coefficient  (see 
Eq.  B.l) . * 
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We  have  pointed  out  that  the  symmetry  of  the  source 
determines  the  partitioning  of  energy  into  the  various 
terms  of  the  series.  But  for  the  cratering  calculations 
the  displacement  field  is  computed  on  only  one  hemisphere. 

It  is  necessary  to  assume  values  on  the  other  hemisphere  to 
compute  the  multipole  coefficients.  The  way  this  assumption 
is  made  can  have  a substantial  influence  on  the  solution. 

Let  us  now  look  more  closely  at  the  problem  of  im- 
mediate interest  as  we  describe  our  modification  of  the 
multipolar  expansion  technique  for  this  application.  In 
Figure  B.2  the  geometry  and  coordinate  system  for  a typical 
calculation  is  shown.  The  explosion  is  at  a depth,  h,  in 
a homogeneous  half space.  The  radius,  R , on  which  the  curl 
and  divergence  of  the  displacement  field  is  monitored  varies 
from  2.3  h to  24.5  h for  the  twelve  calculations  to  be 
studied,  though  it  is  between  5 h and  10  h for  most  of  them. 

Now,  how  are  we  to  select  values  for  the  curl  and 
divergence  in  the  upper  hemisphere  (9e(0,ir))  in  order  to 
make  the  data  compatible  with  our  multipolar  representation? 
The  simplest  and  most  natural  choice  is  to  make  the  data 
in  the  upper  hemisphere  depend  on  that  in  the  lower  hemi- 
sphere in  some  simple  fashion. 

In  Appendix  C we  discuss  the  mathematical  consequence 
of  certain  kinds  of  source  symmetry.  There  are  two  cases 
of  particular  interest;  that  of  axisymmetry  plus  vertical 
symmetry  which  is  characterized  by 


X(4)  (ir-e)  = x(4)  (6)  , 


X(1)  (t-9)  = -X(1)  (0)  , 


x(2)  (IT- 0)  = -X(2)  (9)  , 


and  is  discussed  in  Section  C.2.2  and  that  of  axisymmetry 
and  vertical  antisymmetry: 
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Figure  B.2.  The  geometry  and  coordinate  system  for  the 

cratering  calculations.  The  solution  is  inde- 
pendent of  the  azimuthal  coordinate. 
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X(4)  (tt-0)  = -x(4)  (9)  , 


x(1)  (TT-e)  = x(1)  (Q)  / 


X(3)  (ir-0)  = X{2)  (9)  , 

discussed  in  Section  C.2.3.  The  non-zero  multipole  coeffi 
cients  for  these  two  cases  are  listed  in  Table  C.3.  In 
Figure  B.l  we  showed  the  radiation  patterns  for  each  term 
in  the  multipolar  expansion  of  the  divergence  or  P wave 
portion  of  the  displacement  field.  If  we  assume  vertical 
symmetry  the  even  order  terms  (monopole,  quadrupole,  etc.) 
will  be  non-zero.  If,  on  the  other  hand,  vertical  anti- 
symmetry is  assumed,  the  odd  order  terms  are  non-zero. 

The  radiation  patterns  for  the  lower  order  terms  are 

important  because  these  terms  dominate  the  response  at 

long  periods.  In  fact,  we  can  be  certain  that  the  very 

first  term,  monopole  or  dipole  depending  on  the  symmetry 

chosen,  will  be  the  only  term  that  matters  for  computing 

far-field  surface  waves  and  M . Therefore,  let  us  look 

s 

at  the  radiation  patterns  for  the  S wave  coefficients. 

From  Table  C.3  and  Eq.  (A. 4)  we  see  that  the  radiation  pat 
terns  are  given  by  the  Legendre  functions  P^(cos0)  with  i 
even  for  vertical  symmetry  and  i odd  for  vertical  anti- 
symmetry. These  patterns  are  plotted  in  Figure  B.3. 

How  are  we  to  choose  the  most  appropriate  vertical 
symmetry?  The  first  suggestion  that  comes  to  mind  is  to 
examine  the  angular  variation  of  the  computed  data  at 
various  time  points.  Typical  radiation  pattern  plots  of 
the  data  are  shown  in  Appendix  D where  we  describe  several 
calculations  in  considerable  detail.  We  see  that  these 
radiation  pattern  plots  are  not  much  help.  There  are  at 
least  four  reasons  why  this  is  so: 
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Figure  B.3.  The  radiation  patterns  for  the  first  several 
terms  of  the  curl  for  an  axisymmetric  source. 
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1.  The  calculated  divergence  and  curl  include  the 
influence  of  Rayleigh  waves  at  late  times. 

We  want  to  choose  a symmetry  that  minimizes 
the  effect  of  this  contribution. 

2.  The  numerical  calculations  are  noisy.  Our  multi- 
polar source  representation  should  act  as  a 
smoothing  filter  that  pulls  the  important  signal 
from  this  noisy  data. 

3.  The  calculations  are  really  not  yet  complete; 
that  is,  energy  is  still  propagating  through  the 
elastic  radius  at  the  last  time  point.  We  are 
mainly  interested  in  imposing  the  right  symmetry 
at  long  times  as  the  high  frequency  solution  is 
much  less  dependent  on  this  choice. 

4 . Our  real  interest  is  in  the  radiation  patterns  for 
the  far-field  component  of  the  source  (that  portion 
decaying  as  R ^) . The  patterns  for  the  computed 
data  are  strongly  influenced  by  near-field  terms. 

We  can  attempt  to  gain  some  insight  with  analytical 
methods.  Among  the  class  of  simple  problems,  perhaps  the 
closest  to  the  creating  calculations,  at  least  at  long  times, 
is  the  case  of  a vertical  point  force  on  the  surface  of  a 
homogeneous  half space.  The  far-field  radiation  patterns 
for  the  divergence  and  curl  are  shown  for  this  source  in 
Figure  B.4.  We  see  that  the  P wave  pattern  is  very  much 
like  that  of  a dipole  (see  Figure  B.l).  This  suggests  that 
vertical  antisymmetry  should  work  very  well  for  the  P wave 
portion  of  the  radiation  field.  For  the  S waves  vertical 
symmetry  would  probably  be  a better  choice  since  the  curl 
radiation  pattern  in  Figure  B.4  looks  more  like  a quadrupole. 
But  from  Appendix  C,  especially  Tables  C.l  and  C.2,  we  see 
that  choice  of  different  symmetries  for  the  curl  and  diver- 
gence leads  to  inconsistencies  in  the  displacement  field. 
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Figure  B.4.  Far-field  radiation  patterns  for  the  divergence 
and  curl  due  to  a point's  load  on  the  surface 
of  a homogeneous  half space. 
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Therefore,  we  choose  vertical  antisymmetry  and  represent 
the  field  as  a sum  of  terms  of  order  1,  3,  5,  .... 

It  turns  out  that  we  need  not  be  overly  concerned  with 
the  accuracy  of  our  S wave  representation  anyway.  This  is 
because  the  S waves  must  be  many  times  larger  than  the  P 
waves,  a most  unlikely  circumstance,  to  have  an  appreciable 
effect  on  either  m^  or  Ms-  This  point  is  addressed  in  the 
main  body  of  the  report.  Sections  VI  and  VII. 

What  about  the  contribution  of  locally  generated  Rayleigh 
waves  to  our  solution?  First,  we  point  out  that  the  dis- 
placements due  to  the  Rayleigh  waves  will  only  affect  the 
solution  in  the  region  near  the  free  surface.  Referring 
to  Figure  B.2,  we  don't  expect  the  Rayleigh  waves  to  be  of 
importance  for  takeoff  angles,  t,  less  than  say,  60°. 

Second,  the  Rayleigh  wave  has  a much  greater  effect  on  the 
S wave  or  curl  portion  of  the  field  than  on  the  P wave 
portion.  These  points  were  validated  by  an  elastic  test 
problem  where  we  applied  our  procedures  to  the  field  from  a 
point  source  on  the  surface  of  a homogeneous  half space.  Once 
again,  we  expect  the  errors  to  be  greatest  in  the  S wave 
terms  and  at  long  periods.  Fortunately,  the  P wave  terms 
are  dominant  for  m,  and  M determinations. 
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APPENDIX  C 

PROPERTIES  OF  AN  AXIALLY  SYMMETRIC  SEISMIC  SOURCE 

C.l  SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  AND  POTENTIAL 

FIELDS  FROM  AN  AXIALLY  SYMMETRIC  SOURCE 

In  studies  of  the  ground  motion  from  many  sources  the 
calculation  is  greatly  simplified  if  the  seismic  source 
possesses  the  axial  symmetry.  Teleseismic  ground  motion 
from  many  complicated  seismic  sources  can  then  be  computed 
by  merging  a 2D  nonlinear  finite  difference  source  calcula- 
tion with  the  elastic  wave  propagation  methods  of  theoreti- 

3 

cal  seismology.  This  has  been  accomplished  routinely  at  S 
by  means  of  an  equivalent  elastic  source  representation  of 
the  seismic  source.  The  theoretical  foundation  for  this 
representation  appears  in  a section  of  a recent  paper  by 
Bache  and  Harkrider  [1976] . This  section  is  reproduced  as 
Appendix  A of  this  report. 

Many  sources  not  only  possess  axial  symmetry  but  also 
exhibit  an  additional  physical  symmetry  with  respect  to  the 
plane  normal  to  the  axis  of  symmetry.  We  will  refer  to 
two  special  cases  of  this  symmetry,  that  of  vertical  sym- 
metry and  that  of  vertical  antisymmetry.  Then  for  axially 
symmetric  sources,  three  types  of  symmetry  will  be  dis- 
cussed : 

1.  Axial  symmetry. 

2.  Axial  symmetry  and  reflection  symmetry. 

3.  Axial  symmetry  and  reflection  antisymmetry. 

A seismic  source  exhibits  axial  symmetry  if  the  source 
is  invariant  under  rotation  about  a given  axis,  say  the 
z-axis.  As  a consequence  of  this  symmetry,  the  particle 
displacement  field  must  be  invariant  with  respect  to 
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rotations  about  the  z-axis  and  consequently  have  no  azimuth 
component.  Since  the  field  is  independent  of  the  azimuthal 
coordinate,  <t>,  the  problem  can  be  reduced  to  two  dimensions. 
This  results  in  a substantial  reduction  in  the  computational 
effort.  In  general,  for  an  axially  symmetric  source  the  2D 
finite  difference  calculation  can  be  limited  to  an  area  en- 
closed by  a semi-circle  (0°  <_  0 <_  180°).  In  the  special 
case  where  the  source  exhibits  either  vertical  symmetry  or 
antisymmetry,  the  2D  source  calculation  can  be  further 
limited  to  a single  quadrant  (0°  <_  0 <_  90°).  The  calcula- 
tion of  the  multipole  coefficients  representing  the  equiva- 
lent elastic  source  requires  a double  integration  in  the 
spherical  coordinates  <t>  and  9 over  a spherical  surface  in 
the  elastic  region  that  completely  surrounds  the  nonlinear 
zone  containing  the  seismic  source.  For  an  axially  symmetric 
source,  however,  this  surface  integration  can  be  reduced  to 
a single  integral  along  a circular  arc  in  the  elastic  region. 
This  linkage  between  the  2D  finite  difference  source  calcula- 
tion and  the  elastic  wave  propagation  methods  of  theoretical 
seismology  is  accomplished  by  the  MULTEES  series  of  computer 
programs  which  calculate  the  multipole  coefficients  for  the 
equivalent  elastic  source.  The  2D  source  calculation  must 
be  carried  out  into  the  elastic  region  where  time  histories 
can  be  saved  for  the  displacement  potentials,  which  are  de- 
fined below  in  terms  of  the  divergence  and  curl  of  the  dis- 
placement field.  These  time  histories  are  generated  at  a 
set  of  monitoring  stations  that  are  either  (1)  on  the  arc 
of  fixed  radius  in  the  elastic  region  for  a 2D  calculation 
in  polar  coordinates  (r,0);  or  (2)  in  the  neighborhood  of 
an  arc  in  the  elastic  region  for  a 2D  calculation  in  rec- 
tangular coordinates  (Z,Y).  MULTEES  performs  all  of  the 
operations  required  to  transform  the  2D  variables  for  the 
calculation  of  the  multipole  coefficients,  in  the  latter  case 
including  a 2D  spatial  interpolation  to  evaluate  saved 
variables  on  a circular  arc.  The  MULTEES  programs  are 
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described  in  Bache,  et  al . [1975].  A typical  2D  section  of 

an  axially  symmetric  source  at  fixed  azimuthal  angle  <J>  is 
illustrated  in  Figure  C.l.  The  Z-axis  is  the  axis  of  sym- 
metry and  the  X-Y  plane  is  a possible  plane  of  symmetry. 

The  relationship  between  Cartesian  (1,2,3),  spherical  (R,9,<J>) 
and  2D  finite  difference  code  (Y,Z)  coordinate  system  is  also 
illustrated. 

The  symmetry  properties  of  the  particle  displacement 
field  will  be  used  to  derive  the  symmetry  properties  of 
the  scalar  (dilatational ) displacement  potential  and  the 
vector  (rotational)  displacement  potential.  The  vector  and 
scalar  potentials  are  defined  in  terms  of  the  displacement 
field  u(r,t)  by, 

X = i V x u , 

(C.l) 

X(4)  = V • u. 

The  symmetry  properties  of  the  displacement  potentials  will 
then  be  applied  to  derive  the  nonzero  multipole  coefficients 
for  the  equivalent  elastic  source  representation  of  axially 
symmetric  sources  in  each  of  the  three  classes  defined 
above . 

For  an  axially  symmetric  source  the  displacement  field 
is  invariant  to  rotations  about  the  Z-axis  and  has  no  azi- 
muthal component.  As  a consequence,  the  dilatational  poten- 
(4 ) 

tial  x must  be  independent  of  <f>  and  the  rotational 
potential  x must  have  only  a single  non-vanishing  component 
in  spherical  coordinates,  x^/  which  is  also  independent  of 
4> . For  an  axially  symmetric  source  the  non-vanishing  dis- 
placement and  potential  components  in  spherical  coordinates 
are 
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Figure  C.l.  Coordinate  system  for  an  axially  symmetric 
source  at  fixed  azimuthal  angle.  The  inset 
illustrates  the  relation  between  the  Cartesian 
(1,2,3),  spherical  (R,9,<j>)  and  2D  finite  dif- 
ference code  (Y , Z ) coordinate  systems. 
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Ur(R,9,<j>,t)  = ur  (R,6,t) 

UQ  (R,0,<!>,t)  = u (R, 9 , t)  (C  . 2 ) 

(R,d,<p,t)  = 0 
xr  (R/9»4>/t)  = 0 

Xe  (R,e,<p,t)  = 0 (C . 3 ) 

X^(r,9f<J>#t)  = (R,  9 , t ) 

The  Cartesian  coordinates  of  the  displacement  potentials  are 
X ^ (R,  9,<j)  ,t)  = - sinc(>  x<+)(Rf0»t) 

(2) 

X ; (R,  9 , $ , t)  = cos<}>  x^  (R,  6 , t) 

(C . 4 ) 

X (3)  (R/  9 , , t ) = 0 
( 4 ) 

x (R,e,4),t)  = x4  (R,e,t) 


Note  that  if  the  2D  finite  difference  source  calculation  is 
carried  out  in  the  YZ  plane,  then  one  can  identify  the  azi- 
muthal component  x^  as, 

(0  , Y,  Z , t)  J tcurl  u]  = “ J ("9Y  “ Tz)  ' 

Because  of  the  reduction  in  the  number  of  non-vanish- 
ing spherical  components  of  u and  x*  the  symmetry  properties 
that  are  characteristic  of  an  axially  symmetric  source  can 
be  first  expressed  for  the  spherical  components  and  then  de- 
rived for  Cartesian  coordinates  by  means  of  the  coordinate 
transformation  for  any  vector  A, 
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sin0  cosiji  A + cos9  cosd>  A„ 
r 0 

sin0  sin<}>  A^  + cos0  sincf>  Ag 


- sin<f)  A , 


+ COS(J)  A 

<P 


(C . 5 ) 


A,  = COS0  A 
3 r 


sin0  A0 . 


The  symmetry  properties  of  an  axially  symmetric  source 
are  given  for  the  particle  displacement  field  in  Table  C.l. 
The  symbols  u,  v,  and  w represent  the  X,  Y and  Z components 
of  the  displacement  field  u.  Values  for  these  components 
in  the  second,  third  and  fourth  octants  are  related  to  their 
corresponding  value  in  the  first  octant.  If  either  case  of 
vertical  symmetry  obtains,  an  additional  relation  exists  for 
points  below  the  Z = 0 plane  (fifth  octant) . Consequently, 
displacement,  velocity  and  potential  components  at  all  points 
on  the  surface  of  a sphere  in  the  elastic  region  can  be  de- 
rived from  the  time  histories  saved  on  the  circular  arc 
0°  < 0 < 180°  for  axial  symmetry  or  on  the  arc  0°  <_  0 <_  90° 
for  axial  plus  vertical  (anti)  symmetry. 


The  symmetry  properties  of  the  dilatational  and  rota- 
tional displacement  potential  can  be  derived  from  those  for 
the  displacement  field.  The  results  of  this  derivation  are 
given  in  Table  C.2.  One  example  illustrates  the  method 
employed.  From  the  definition  of  the  rotational  potential, 
Eq.  (C.l)  and  the  symmetry  properties  for  vertical  symmetry, 


XU)  = j [curl  u] 


1_  / _ 3v  \ 
2 \3y  " 3z ) 


(C.  6) 


X (1) (x,y,-z) 


1 /3 ( — w ) 

2 \ 3y 


v 

(-Z) 


(x,y , z) 


As  a consequence  of  the  symmetry  properties,  the 
following  boundary  conditions  will  apply: 
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TABLE  C.l 

SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  FIELD  FOR  AN  AXIALLY 

SYMMETRIC  SOURCE 


Field  i 

components  in 

the  2nd,  3rd, 

4th  and  5th 

octants  (column 

2)  are 

related  to  their  values 

in 

the  1st  octant  (columns  3- 

5)  . 

(u,  V 

,w  = X,  Y and 

Z components 

of  particle 

displacement) . 

Axial 

Field 

Symmetry 

Vertical 

Vertical 

Octant 

Component 

Only 

— 

Symmetry 

Antisymmetry 

5 

ur  (r,n-0) 

★ 

+ur  (r , 0) 

-ur  (r,0) 

u g (r , iT-0 ) 

★ 

1 

c 

CD 

CD 

+uQ  (r , 0 ) 

5 

u (x,y ,-z) 

* 

+u (x,y , z) 

-u (x,y,z) 

v (x,y,-z) 

* 

+v (x,y, z) 

-v  (x,y , z) 

w(x,y,-z) 

★ 

-w  ( ) 

+w(  ) 

2 

u (-x,y , z) 

-u (x,y. 

z) 

-u (x,y ,z) 

-u  ( ) 

v (-x ,y , z ) 

+v( 

) 

+v(  ) 

+v(  ) 

w (-x,y, z) 

+w  ( 

) 

+w  ( ) 

+w  ( ) 

3 

u (-x,-y ,z) 

-U( 

) 

-u(  ) 

-U  ( ) 

v (-x,-y, z) 

-V( 

) 

-V  ( ) 

-V(  ) 

w (-x , -y , z ) 

+w  ( 

) 

+w  ( ) 

+w  ( ) 

4 

u(x,-y,z) 

+u( 

) 

+u(  ) 

+u(  ) 

v (x,-y, z) 

-V  ( 

) 

-V(  ) 

-v(  ) 

w (x, -y , z) 

+w  ( 

) 

+w  ( ) 

+w  ( ) 

★ 

Relationship  undefined 


88 


R-3119 


TABLE  C.2 

SYMMETRY  PROPERTIES  OF  THE  DISPLACEMENT  POTENTIALS  FOR  AN 

AXIALLY  SYMMETRIC  SOURCE 

Field  components  in  the  2nd,  3rd,  4th  and  5th  octants  (column 
2)  are  related  to  their  values  in  the  1st  octant  (columns  3- 
5)  . 

(X(1),  X(2),  X^  = X,  Y and  Z components  of  rotation  dis- 
placement potentials) . 


Axial 

Field 

Synroetry 

Vertical 

Vertical 

Octant 

Component 

Only 

Syirrietry 

Antisymrietry 

5 

(r,Tr-0) 

★ 

-yr'e> 

+X^(r,0) 

5 

X(1) (x,y,-z) 

★ 

-x(1)  (x,y,z) 

+X(1)(x,y,z) 

X(2)  (x,y,-z) 

* 

-x(2,(  ) 

+x <2)  < ) 

x(4)(  ) 

* 

+X(4) (x,y,z) 

-x<4) ( ) 

2 

X(1)  (~x,y,z) 

+X(1)  (X,y,: 

z) 

+X(1> (x,y,z) 

+Xa>  (x,y,z) 

X(2)  (~x,y,z) 

-x(2>( 

) 

+x(2,(  ) 

V2)<  ) 

x(4)(  ) 

V4,< 

) 

y4)<  ) 

+x(4)<  ) 

3 

X(1)  (~x,-y,z) 

-x(1)( 

) 

-x(1)<  ) 

-x(1><  ) 

X(2) (-x,-y,z) 

-x<2)( 

) 

-x(2)(  ) 

-x(2)(  ) 

x(4)(  ) 

+x(4)< 

) 

V4,<  ) 

+x(4,(  ) 

4 

X(1)  (x,-y,z) 

-x(1)( 

) 

-x(1)(  ) 

-x(1)(  ) 

X(2)  (x,-y,z) 

V2)< 

) 

V2><  ) 

+x<2)  ( ) 

x(4)(  ) 

♦x(4)  ( 

) 

+x!4,<  ) 

+x<4)  ( ) 

Relationship  undefined 
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1 • Axial  Symmetry  Only 

On  the  Z-axis  corresponding  to  x = y - 0 or 

0 = 0,  x = 0 or 

X(1)  (R,0,t)  = 0, 

X(2) (R,0,t)  = 0.  (C . 7 ) 

2 . Axial  Plus  Vertical  Symmetry 

On  the  XY  plane  (the  plane  of  symmetry)  corres- 
ponding to  z = 0 or  0=90°,  x = 0,  or 

X (1)  (R,ir/2,t)  = 0, 

X(2)  (R,Tr/2,t)  = 0.  (C . 8 ) 

3 . Axial  Plus  Vertical  Antisymmetry 

On  the  Xv  plane  corresponding  to  z = 0 , or 

6 = 90°, 

X(4)  (R,ir/2,t)  = 0.  (C.9) 

C. 2 MULTIPOLE  COEFFICIENTS  FOR  AN  AXIALLY  SYMMETRIC  SOURCE 

The  equivalent  elastic  source  representation  for  any 
arbitrary  seismic  source  is  defined  by  a multipole  expansion 
in  spherical  harmonics.  Calculation  of  the  multipole  coef- 
ficients requires  a double  integration  over  the  surface  of 
a sphere  in  the  elastic  region.  The  symmetry  properties  de- 
fined above  can  be  applied  in  order  to  derive  the  non-vanish- 
ing multipole  coefficients  for  any  axially  symmetric  source. 
While  the  nonzero  coefficients  derived  in  this  manner  define 
the  equivalent  elastic  source  for  a specific  configuration 
in  the  source  coordinate  system,  the  coefficients  for  any 
other  configuration  can  be  obtained  from  this  set  by 
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. 

f 


performing  the  appropriate  transformations  to  the  source 
coordinate  system. 

The  multipole  coefficients  in  the  time  domain  (see 
Bache,  et  al . [May,  1975])  can  be  expressed  as 


jAta  (R't)) 

2tt  it 

* Cta  f f 

0 0 

/cos  m<f>\ 

x P™  (cos6)  sin9d0  J 


'sin  m<t> ' 


d<f> , 


(C.10) 


where 

e (2S.+1)  U-m)  ! 

r = -£ , 

£m  4 tt  (Jl+m)  ! 

with  £q  = 1 and  = 2,  m / 0 and  P™  are  associated  Legendre 
functions;  the  double  integration  is  carried  out  on  a sphere 
of  radius  R in  the  elastic  region  and  a = 1,2, 3, 4 corres- 
ponds to  the  dilatational  (a=4)  and  rotational  potential 
(a=l ,2,3)  . 


C.2.1  Axial  Symmetry  Only 

The  case  in  which  the  seismic  source  exhibits  only 
axial  symmetry  will  be  considered  first.  The  integration 
in  can  be  carried  out  explicitly  by  utilizing  the  follow- 
ing integral  properties; 


(C.ll) 
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2tt 


2tt 


/ 


sinij) 


COSlfl 


cos 

m<{>) 

1° 

dcp  = 

) 

t 

cos 

m4>) 

( * 

m = 1 

cos 

m<J> ) 

( TT 

m = 1 

d<p  = 

I 

\ 

sin 

m<f>) 

( 0 

(C • 12 ) 


(C . 13 ) 


(4) 

Since  the  dilatational  potential  x is  independent  of  <p , 
the  only  nonvanishing  P wave  multipole  coefficient  must  have 
m = 0 , 


Ato>(R't) 


= 2Tr  C 


£0 


/ 


(4) 


(R,9,t)  P^(cos6)  sin0  d9, 


U = 0,1,2,3,...) 


B*m  ,R'U  - 0 


(C . 14 ) 
(C . 15 ) 


Using  the  relationship  (Eq.  (C.4))  between  the  two  nonzero 
rotation  potential  components  and  x^  and  the  azimuthal 

potential  the  rotation  potentials  can  be  reduced  to  the 

following  single  integrals. 


Bta  (E't» 


2 TT 


= Clmf  x (R,Q,t)  P™  (cos0)  sin0  d9 


2ir 


cos  m<J> 


(-sin<f>)  'd<j>, 

n ( sin  m<J> ) 


(C . 16 ) 


(R.t)  « 0 


(C . 17 ) 
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IT 

(P,t)  = - tt  CJll  J X({)(R»e/t)  (cos6)  sin6  d6, 

0 

(4  = 1,2,3, .. .)  (C . 18 ) 

IT 

C4m  f X*(R,e,t)  P™  (cose)  sine  d6 
0 


jAim>(R't)| 

I B«>  (R.t)l 


/<  cos  m<fr  j 

cos#  i d<(>,  (C.19) 

>sin  m<}> ) 


ir 

A^}  (R,t)  = C£1  J X(J)(R,0,t)  pj  (cose)  sine  d6 , 

0 


U = 1,2,3,...) 

(C.20) 

Bta  <R't)  * 0 

(C. 21) 

> 

3 “ 

50 

ft 

II 

o 

(C.22) 

<*'«  * 0 

(C. 23) 

This  completes  the  derivation  of  the  nonvanishing  multipole 
coefficients  for  an  axially  symmetric  source  that  does  not 
exhibit  any  other  special  symmetry.  One  result  of  this  re- 
duction is  that  only  a single  S wave  rotational  potential 
needs  to  be  calculated,  since 

(R,t)  = - B^^  (R,  t)  . (C . 24 ) 
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The  divergence  and  curl  of  the  displacement  field  are 
monitored  in  the  2D  finite  difference  source  calculation 
at  stations  on  the  circular  arc  of  radius  R in  the  elastic 
region.  If  the  coordinates  used  in  the  2D  source  calcula- 
tion are  (Y,Z)  for  a rectangular  grid  or  (R,9)  for  a polar 
grid  (see  Figure  C.l),  one  can  make  the  following  identifi- 
cation, 

X^(R,eft)  = - J CURLD  (R,e,t)  = - J CURLD  (Y,Z,t), 

X(4)(R,9,t)  = DIVD  (R, 0 , t)  = DIVD  (Y , Z , t)  (C.25) 


where 


R = (Z2  + Y2)1/2, 


6 = tan-1  (Y/Z) , 


and  DIVD  and  CURLD  are  the  divergence  and  curl  of  the  dis- 
placement field,  respectively,  as  generated  by  the  source 

(4 ) 

calculation.  By  applying  these  relations  for  x in  Eq. 

(C.14)  and  for  in  Eq.  (C.19),  the  numerical  integration 
<P 

in  9 is  carried  out  in  a straightforward  fashion. 


C.2.2  Axial  and  Vertical  Symmetry 

Whenever  the  seismic  source  exhibits  symmetry  or  anti- 
symmetry with  respect  to  reflections  in  the  XY  plane,  the 
integration  in  0 can  be  reduced  to 

ir/2 

f [ 1 + <-)4+m]  x(C°  <R,9,*,t)  P ” (cos0)  sin0  d0 , 

0 (C. 26) 

Since  P™  (cos(it-0))  * (-)  ^+m  p^1  (cos0)  ; the  positive  (nega- 
tive) sign  in  the  bracket  corresponds  to  the  case  where 
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( cx ) 

X is  symmetric  (antisymmetric)  with  respect  to  reflection 
in  the  XY  plane.  For  the  case  of  axial  and  vertical  sym- 
metry, 


X(4)  (ir-0)  = + X(4) 


(9)  , 


X(1)  Or-6)  = - xU)  (9)  , 


X(2)  (*-9>  = - X(2)  (9)  . 


(C. 27) 


Consequently,  the  nonvanishing  P wave  multipole  coefficients 
are  reduced  to 


tt/2 

A^4)  (R,t)  = 4 7r  ClQ  J x (4)  (R/  9 /t)  P£  (cos  0 ) sine  d6, 

0 

(2.  = 0,2,4.  . .)  (C . 28) 

Note  that  the  integration  in  <p  first  eliminates  all  terms 
except  those  with  m = 0;  the  9 integration  eliminates  all 
odd-order  terms.  The  nonvanishing  S wave  multipole  coef- 
ficients are  reduced  to 


tt/2 

(R,t)  - - 2ir  C£1  J X(J)(R/9,t)  pj  (cos9)  sin6  d0 , 

0 

U = 2,4,6,...)  (C. 29) 

ir/2 

AJll^  (R't)  ~ 2lT  J X^(R#9ft)  P^  (cos0)  sin9  d6, 

0 

(*■  * 2,4,6,...)  (C . 3 0 ) 
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In  summary,  the  calculation  of  the  multipole  coeffi- 
cients for  a source  having  axial  as  well  as  vertical  symmetry 

is  greatly  simplified  to  a set  of  integrals  along  an  arc 

(0°  <_  0 <_  90°)  at  a radius  R in  the  elastic  region.  Only  a 

single  term  is  found  to  be  nonzero  for  each  even  order 

multipole,  namely  that  corresponding  to  m = 0 for  P wave 
terms  and  m = 1 for  S wave  terms.  The  leading  term  in  the 

equivalent  elastic  source  representation  corresponds  to  a 

(4 ) 

P wave  monopole  (£=0)  AqQ  , followed  by  both  P wave  and  S 
wave  quadrupole  terms  (£=2) , A^j^ , > A2i^  * Vertical 

symmetry  is  often  referred  to  as  "quadrupole"  symmetry. 
Multipole  coefficients  for  complex  axially  symmetric  sources 
were  first  calculated  and  applied  in  teleseismic  studies  by 
Cherry,  et  al.  [May,  1975]  . 

C.2.3  Axial  and  Vertical  Antisymmetry 

For  an  axially  symmetric  source,  which  also  exhibits 
vertical  antisymmetry  with  respect  to  reflections  in  the  XY 
plane,  the  relations 

X(4)  (tt-9)  = - X(4)  O)  , 

X(1)  (w-9)  = + X(1)  (9)  , (C . 31) 

X(2)  (ir-9)  = + X(2)  (9)  , 

Can  be  employed  in  Eq.  (C.25)  to  derive  the  following  non- 
vanishing multipole  coefficients, 

ir/2 

A^4)  (R,t)  * 4 ir  C£Q  J x(4)(R,9,t)  (cos0)  sine  d6 , 

0 

(l  = 1,3,5,...)  (C. 32) 
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tt/2 

B^;  (R,t)  = - 2tt  CZ1  J X^(R»9ft)  pJ  (cos9)  sin0  d9, 

0 

(Z  - 1,3,5, ..  .)  (C. 33) 

A^}  (R, t)  = - B^}  (R, t ) . 

(Z  = 1,3,5,...)  (C . 34 ) 

The  leading  terms  in  the  equivalent  elastic  source  representa- 

(4 ) 

tion  corresponds  to  a P wave  dipole  (£=1)  A'  and  an  S wave 
(1)  (2) 

dipole  U-l)  and  A^  . Vertical  antisymmetry  is  often 

referred  to  as  "dipole"  symmetry.  There  is  no  monopole  or 
quadrupole  terms  in  such  a case.  The  nonzero  multipole 
coefficients  are  summarized  in  Table  C.3. 


C. 3 FAR- FI ELD  DISPLACEMENT  SPECTRA  FOR  AN  AXIALLY  SYMMETRIC 

SOURCE 

For  the  determination  of  the  teleseismic  signature  for 
an  axially  symmetric  source,  only  the  far-field  portion  of 
the  displacement  field,  which  includes  only  those  terms  de- 
caying with  R ^ in  the  multipole  expansion,  is  required. 

The  displacement  spectrum  can  be  expressed  in  terms  of  the 
potential  spectrum  by 


u (R,u>) 


v x(RfW) 


* X (R/W) r 


where 


(C. 35) 


X(4)(R,w)  (w)  P£  (cos 6 ) h^2)  (kpR),  (C.36) 

Z=0 
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(R,w) 


y]  B<2'  (u)  pj  (cos6)  Sin4>  hjJ2)  (kgR)  , 

£=1 

(C. 37) 


X(2)(Rr^)  = ^2  Ajj2)  ((d)  pj  (COS0)  COS<|>  hj[2) 
£■1 


(ksR)  , 

(C- 38) 


X(3) (RrU)  = 0, 


with 


, -(1) 
- cos^)  x 


k = w/a,  k = o)/6 , where  a and  6 are  the  P wave  and  S wave 
P ® 

body  wave  velocities,  respectively;  h^  are  spherical  Hankel 
functions  of  the  second  kind;  and  the  vector  R has  as 
spherical  components  R , © , <j> . The  multipole  coefficients  are 
calculated  by  taking  the  Fourier  transform  as  defined  by 


' h~«y/k  dt.  (C.39) 

n 'Ka  -00 

The  far- field  spherical  components  of  the  displacement  spec- 
trum are 
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u (R,co) 

<t>  ~ 


(C. 42) 


where  in  the  far-field  and  radial  displacement  results  only 
from  P wave  terms  involving  and  the  transverse  displace- 

ment is  only  excited  by  S wave  terms  involving  x- 

In  general,  for  an  axially  symmetric  source  the  far- 
field  displacement  spectrum  reduces  to 
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" f the  source  exhibits  axial  plus  vertical  symmetry,  the  lead- 
ing terms  for  the  displacement  spectrum  are 


uR  ( R , co ) 


FF 


k R 
P 


-ikpR  ( 

(co)  - (co)  P2  (cose) 


+ A^  (co)  P4  (cose) 


+ • • • 


(C. 45) 


ue(?'u)  FF 


”iksR  ( ?(1)  ,..x  „1 


k R 
s 


(<»>)  p2  (cose) 


+ B^  (co)  pj  (cose) 


+ * * * I * 


(C. 46) 


100 


R-3119 


For  the  case  of  axial  symmetry  plus  vertical  antisymmetry, 
the  leading  terms  for  the  displacement  spectrum  are 
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APPENDIX  D 

DESCRIPTION  OF  A TYPICAL  CALCULATION 


Teleseismic  ground  motion  is  determined  by  linking  the 
source  calculations  performed  by  Applied  Theory,  Inc.  (ATI) 
with  elastic  wave  propagation  methods.  This  has  been  ac- 
complished numerically  by  means  of  the  MULTEES  series  of 
general  purpose  computer  routines  developed  at  S3.  MULTEES 
calculates  the  MULTipole  coefficients  for  the  Equivalent 
Elastic  Source  representation  for  any  arbitrary  numerical 
source,  which  has  been  carried  into  the  small  displacement 
elastic  region.  MULTEES  accepts  the  output  generated  in 
either  a rectangular  or  spherical  (polar)  coordinate  grid 
system  by  any  two  or  three  dimensional  finite  difference 
source  code. 

A brief  outline  of  the  key  steps  in  the  computational 
procedure  will  be  presented  in  order  to  illustrate  selected 
features  of  a typical  calculation.  Three  numerical  proce- 
dures will  be  described,  namely: 

1.  Preprocessing  of  the  ground  motion  fields  and 
generation  of  the  displacement  potential  fields 
(MULTEES. M2) ; 

2.  Calculation  of  the  multipole  coefficients  in  the 
time  domain  for  the  equivalent  elastic  source 
(MULTEES. M4) ; 

3.  Calculation  of  the  Fourier  transform  of  the 
multipole  coefficients  (MULTEES .M5) . 

The  corresponding  computer  program  applied  in  each  step  is 
indicated  in  parenthesis.  Once  the  set  of  multipole  coef- 
ficients has  been  generated  at  selected  frequencies,  far- 
field  displacement  spectra  and  ground  motion  can  be  cal- 
culated. 
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D. 1 PREPROCESSING  OF  THE  2D  SOURCE  OUTPUT 

D.1.1  Raw  2D  Source  Calculation  Output 

The  configuration  for  a near  surface  explosion  source 
with  a depth  of  burial  given  by  DOB  is  illustrated  in  Figure 
D.l.  For  each  of  the  sources  calculated  by  ATI,  the  2D  out- 
put data  for  the  ground  motion  fields  consist  of  the  follow- 
ing quantities  at  190  stations  for  each  time  cycle: 

T Time  (sec) 

UH  Horizontal  component  of  velocity  (cm/sec) 

UV  Vertical  component  of  velocity  (cm/sec) 

SIGH  Horizontal  component  of  stress  (dynes/cm^) 

SIGV  Vertical  component  of  stress  (dynes/cm^) 

SIGPHI  Azimuthal  component  of  stress  (dynes/cm2) 

SIGHV  In-plane  component  of  shear  stress  (dynes/cm^) 
DIVD  Divergence  of  the  displacement  field 
CURLD  Curl  of  the  displacement  field. 

The  direction  of  the  curl  is  that  of  the  cross  product 

A A A A 

u^  * u^,,  where  u^  and  uv  define,  respectively,  a positive 
horizontal  direction  and  the  upward  vertical  direction.  The 
first  91  monitoring  stations  were  located  at  1-degree  inter- 
vals, starting  at  the  horizontal  axis,  on  a circle  of  radius 
R^ , whose  center  lies  at  the  ground  surface  on  the  axis  of 
symmetry.  The  second  group  of  91  stations  were  also  at  1- 
degree  intervals  but  at  a radius  R2*  There  are  also  8 sta- 
tions at  intervals  of  2.5  degrees  at  a radius  of  R^ , start- 
ing at  the  axis  of  symmetry,  where  R^  > R2  > • The  relation- 

ship between  the  spherical  coordinates  (R,6)  and  the  rectangu- 
lar coordinates  (Z,Y)  used  in  the  2D  finite  difference  cal- 
culation is  also  illustrated  in  Figure  D.l.  The  Z and  Y axes 
represent  the  vertical  and  horizontal  directions,  respec- 
tively. 
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The  symmetry  properties  of  the  displacement  and  po- 
tential fields  of  an  axially  symmetric  source  are  discussed 
in  Appendix  C.  The  derivation  of  the  non-zero  multipole 
coefficients  for  an  axially  symmetric  source,  which  is 
also  assumed  to  have  vertical  (reflection)  symmetry  or 
antisymmetry,  is  given  in  Appendix  C.2.  There  it  is  shown 
that  the  calculation  of  the  multipole  coefficients  in  the 
time  domain  can  be  reduced  to  a numerical  integration  of 
these  potentials  over  a single  quadrant  of  the  circle  of 
radius  R,  namely  from  9 = 90°  to  0 = 180°.  The  required 
Cartesian  displacement  potentials  are  defined  in  terms  of 
the  2D  source  output  variables  (Y,Z)  as  follows: 

X(1)  (R,  9 , <J>=it/2  , t)  = j CURLD  (Y,Z,t) 


(4) 

X (R,  0 , t)  = DIVD  (Y , Z , t ) 


(D.l) 


where 


2 2 
R = (Y  + Z ) 


1/2 


0 = tan-1  (Y/Z), 

and  DIVD  and  CURLD  are  the  divergence  and  curl  output  source 
variables  from  the  2D  finite  difference  calculation. 

The  first  computational  step  in  a typical  calcula- 
tion consists  of  a detailed  examination  of  the  raw  2D  source 
output  data.  This  is  accomplished  by  a series  of  computer 
runs  applying  various  facilities  of  the  MULTEES  program  M2, 
including  such  operations  as 

1.  Plotting  versus  time  all  raw  output  variables 
at  selected  stations; 

2.  Plotting  versus  angle  all  raw  output  variables 
at  selected  times; 
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3.  Listing  the  complete  set  of  raw  output  data  at 
selected  times. 

For  monitoring  purposes,  simple  "printer  plots",  which  are 
generated  routinely  by  the  MULTEES  programs,  are  most  con- 
venient and  appropriate.  In  addition  to  the  eight  raw  out- 
put variables  defined  above,  radial  (R)  and  colatitude  (e) 
components  of  particle  velocity  and  acceleration  are  also 
computed  and  plotted  versus  both  angle  and  time.  The  ob- 
jective of  all  of  this  preprocessing  is 

1.  To  identify  the  general  characteristics  of  the 
raw  source  calculation. 

2.  To  ch^ck  if  the  numerical  solution  has  reached 
or  v-  approaching  a steady- state  condition 

at  the  last  time  cycles. 

.3.  To  identify  any  special  features  or  peculiari- 
ties of  the  numerical  data. 

The  magnitude  as  well  as  the  general  shape  of  the  ground 
motion  fields  were  compared  with  other  2D  axisymmetric  cal- 
culations for  complex  explosion  sources.  As  a result  of 
such  comparisons,  the  features  which  are  unique  to  cratering 
calculations  were  further  highlighted. 

The  angular  variation  of  the  displacement  potentials 
X(1)  and  X(4)  on  the  elastic  radius  at  a few  selected  times 
is  illustrated  in  Figures  D.2-D.4  for  two  typical  cases. 

The  calculations  are  identified  by  number  in  Section  III  of 
this  report.  In  these  graphs  the  horizontal  axis  gives  the 
colatitude  angle  0 in  degrees,  which  measures  the  position 
of  each  of  the  91  monitoring  stations  on  the  elastic  radius 
with  respect  to  the  upward  vertical  position;  the  left  and 
right  boundaries  correspond  to  the  free  surface  (0  = 90°) 
and  the  downward  vertical  direction  (0  = 180°),  respectively. 
Figures  D.2(a),  (b)  and  (c)  display  the  divergence  (with 
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Figure  D.2.  Illustration  of  angular  variation  of  displacement  potential  x ' at 
selected  times  for  Case  1.  Figures  (a),  (b)  and  (c)  display  the  divergence  (with 

overburden  removed)  at  times  of  0.442,  0.471  and  0.506  sec. 


0.508  sec 
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Figure  D.3.  Illustration  of  the  angular  variation  of  the  displacement  potential 
for  Case  2 at  0.508,  1.005  and  1.998  sec. 
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overburden  removed)  for  Case  1 at  times  of  0.442,  0.471 

and  0.506  sec.  These  times  correspond  to  the  arrival  times 

of  the  primary  P wave  signal  at  this  radius  (R^  = 2100  M) 

for  stations  located  at  approximately  6 = 161°,  118°  and 

98°,  respectively.  Figure  D.3  displays  the  dilatational 
(4 ) 

potential  x (divergence  with  overburden  removed)  and 
Figure  D.4  displays  the  rotational  potential  (1/2  curl) 

for  Case  2 at  times  of  0.508,  1.005  and  1.998  sec.  In 
general,  the  angular  variation  of  the  divergence  was  ob- 
served to  be  relatively  smooth  only  for  first  motions,  i.e., 
near  the  primary  P wave  arrival.  At  such  times  the  displace- 
ment and  velocity  patterns  are  also  quite  regular.  However, 
at  most  other  times  the  patterns  for  the  particle  velocity, 
divergence  and  curl  are  noisy  to  a greater  or  lesser  degree. 
The  more  irregular  the  velocity  pattern,  the  noisier  the 
divergence  and  curl  patterns.  In  previous  calculations  per- 
formed at  for  contained  sources,  regular  patterns  are 
observed  throughout  the  time  history;  in  particular,  the  last 
time  cycle  in  such  cases  exhibits  a divergence  and  curl  that 
is  smooth  and  slowly  varying  with  angle. 

D.1.2  Special  Numerical  Features  of  the  Raw  Data 

A detailed  examination  of  the  raw  source  output  data 
identified  two  types  of  special  numerical  features.  The 
first  is  due  to  the  fact  that  the  source  output  variable  DIVD 
included  a contribution  from  the  overburden  pressure.  To 
extract  the  divergence  of  the  displacement  field  from  the 
given  data,  a correction  given  by  pgh/K  must  be  applied, 
where  p is  the  material  density,  g the  gravitational  accelera- 
tion, h the  vertical  distance  of  each  station  below  the  sur- 
face, and  K the  bulk  modulus.  Figure  D.5  illustrates  time 
records  for  Case  12  both  for  the  raw  source  variable  DIVD 
and  the  divergence  after  removing  the  contribution  from  the 
overburden.  About  25  percent  of  the  total  time  history  is 
shown. 
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Figure  D.5.  The  time  history  of  the  divergence  for  the  first  0.6  seconds  of  the  calcula- 
tion for  Case  12.  The  graphs  on  the  left  represent  the  raw  output  (DIVD)  generated  by 
the  source  calculation.  In  the  graphs  on  the  right  the  overburden  correction  has  been 
applied  to  derive  the  divergence  (x*4M  at  0 = 120°  and  0 = 150°. 
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The  second  class  of  discontinuities,  which  were  ob- 
served in  both  the  raw  divergence  (DIVD)  and  curl  (CURLD) 
output  data,  occur  at  times  during  the  source  calculation 
when  rezoning  took  place.  Figures  D.5  and  D.6  illustrate 
this  characteristic  numerical  feature  in  both  the  divergence 
and  curl  time  records  for  Case  12  at  0.351  and  0.562  sec. 
Cases  1,  4,  5,  6,  7,  9,  10,  11  and  12  had  discontinuities 
at  two  different  times,  since  rezoning  generally  occurred 
twice  in  the  course  of  each  source  calculation;  Cases  2,  3 
and  8,  however,  had  discontinuities  at  only  a single  time. 

At  each  discontinuity  there  were  two  divergence  and  curl 
values  given  for  the  same  time  point.  Moreover,  Cases  9 and 
12  were  further  complicated  by  there  being  three  raw  values 
at  the  second  discontinuity  with  one  record  containing  all 
null  values.  This  latter  record  was  simply  skipped.  Ap- 
propriate displacement  potentials  were  extracted  from  the 
raw  data  without  making  any  numerical  adjustment  for  these 
discontinuities . 

D.2  CALCULATION  OF  THE  MULTIPOLE  COEFFICIENT 

Once  the  detailed  examination  of  the  raw  source  data 
and  the  other  required  preliminary  computational  procedures 
described  above  are  completed,  the  MULTEES . M2  program 
generates  the  displacement  potentials  at  91  evenly  spaced 
points  on  the  elastic  radius  from  6 = 90°  to  9 = 180°. 
Complete  time  records  of  the  displacement  potentials  at  a 
single  station  are  illustrated  for  Cases  2 and  12  in  Figures 
D.7  and  D.8.  Two  features  can  be  noticed  in  such  time 
histories.  First,  the  potentials  fail  to  reach  static  values 
at  the  time  the  source  calculation  was  terminated.  Secondly, 
there  appears  to  be  a negative  bias  in  the  divergence  for  the 
sandstone  case  (12).  The  final  time  value  was  2.5  sec  for 
all  cases  except  for  Case  2 where  the  source  time  histories 
were  terminated  at  2.0  sec. 
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Figure  D.7  (Concluded) 
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Figure  D.8  Concluded 
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The  MULTEES  program  M4  uses  the  displacement  poten- 
tials generated  by  M2  to  calculate  the  multipole  coeffi- 
cients in  the  time  domain  according  to  the  symmetry  options 

specified.  The  leading  P wave  and  S wave  terms,  A1  ' (R,t) 

(2)  1U 
and  A^  (R,t),  respectively,  for  axial  plus  vertical  antisym- 
metry (i.e.,  "dipole"  symmetry)  are  displayed  in  Figure  D.7 
and  Figure  D.8.  The  angular  integrations  (see  Equations 
C.32  and  C.33  in  Appendix  C)  are  carried  out  numerically  on 
the  elastic  radius,  which  is  2100  M in  Case  2 and  1200  M in 
Case  12.  A general  feature  observed  for  these  cratering 
calculations  is  the  failure  of  the  multipole  coefficients  to 
reach  a steady-state  value.  Well-behaved  dipole  coefficients 
should  approach  zero  at  late  times.  All  of  the  dipole  coef- 
ficients, however,  in  the  next  processing  step  will  be  set 
equal  to  zero  for  all  times  beyond  the  last  time  point  pre- 
sent in  each  record.  For  higher-order  terms,  the  amplitudes 
are  assumed  to  remain  static  at  the  value  reached  at  the 
last  time  point.  The  length  of  the  time  record  affects  the 
frequency  content  that  can  be  derived  from  the  record. 

The  next  computational  step  uses  the  MULTEES  program 
M5  to  calculate  the  multipole  coefficients  at  a selected 
array  of  frequencies.  One  of  the  first  steps  in  this  cal- 
culation is  to  linearly  interpolate  the  input  time  domain 
multipole  coefficients  to  evenly  spaced  time  points.  This 
process  is  necessary  to  prepare  the  coefficients  for  the  ap- 
plication of  a standard  fast  Fourier  transform  routine. 
Moreoever,  this  removes  a possible  discrepancy  due  to  the 
existence  of  duplicated  records.  For  example,  in  the  raw 
source  output  data  for  Case  2,  a total  of  nine  pairs  of  time 
cycles  were  identified  containing  duplicate  values  for  the 
variables  DIVD  and  CURLD. 

The  fast  Fourier  transform  of  the  multipole  coeffi- 
cient is  taken  for  the  dipole  terms  and  the  transform  of  the 
derivative  of  the  coefficient  is  taken  for  terms  higher  than 
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dipole  for  a total  of  1024  frequencies.  The  frequency  inter- 
val used  in  each  calculation  is  given  as 

Af  = 2048  At  ' 

where 

At  = (T  - T.  )/(n  - 1)  , 
n 1 

T. ,T  = First  and  last  time  record  , 

1 n 

n = Number  of  time  cycles  processed. 

As  an  illustration,  for  Case  2 

Af  = 0.076521  Hz,  T = 2.0015  sec  , 

n 

n = 400,  = 0.15517  sec. 

The  complex  multipole  coefficients  are  transformed  to  ampli- 
tude and  phase.  The  next  step  is  a logarithmic  interpreta- 
tion of  the  amplitudes  and  a linear  interpolation  of  the 
phases  to  extract  values  of  the  multipole  coefficients  for 
a set  of  requested  output  frequencies.  Since  some  of  the 
output  frequencies  requested  may  be  below  the  minimum  Af 

defined  above,  the  M5  program  extrapolates  the  amplitudes 

. -2 
for  the  dipole  *terms  proportional  to  oj  . The  output  from 

this  series  of  calculations  performed  by  MULTEES.M5  is  a 
set  of  multipole  coefficients  for  a specified  array  of  fre- 
quencies. Multipole  coefficients  were  computed  for  the 
twelve  different  cratering  sources  up  to  order  L = 5. 

Figures  D.7  and  D.8(c)  and(f)  display  the  leading  P wave 
and  S wave  terms  A^^  (w)  and  A^^  (w)  , respectively,  for 
Cases  2 and  12.  The  amplitudes  are  plotted  versus  frequency 
on  a log-log  scale  (logarithms  are  in  base  10) . Higher 
order  terms  (2.  = 3 and  2 = 5)  are  illustrated  in  Figures 
D. 9 and  D.10  for  Cases  2 and  12.  It  should  be  pointed  out 
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Figure  D.9  Concluded 
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Figure  D.10  Concluded 
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that  no  damping  or  other  numerical  smoothing  was  applied 
to  the  potential  fields  themselves  nor  to  the  multipole 
coefficients. 

In  summary,  the  final  product  of  the  series  of  compu- 
tational procedures  outlined  above  is  the  set  of  multipole 
coefficients  at  selected  frequencies.  These  define  the 
equivalent  elastic  source  for  each  of  the  cratering  calcula- 
tions. Teleseismic  displacement  spectra  can  then  be  derived 
from  these  coefficients. 
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APPENDIX  E 

THE  SEISMIC  WAVES  DUE  TO  A STRESS  DISTRIBUTION  APPLIED 
AT  THE  SURFACE  OF  A MULTILAYERED  HALFSPACE 


E.l  INTRODUCTION 

We  consider  the  elastic  waves  generated  by  the  appli- 
cation of  an  axisymmetric  stress  distribution  at  the  surface 
of  a medium  consisting  of  plane  elastic  layers.  We  are 
interested  in  both  the  body  and  surface  waves  propagating 
in  the  medium.  Harkrider  [1964,  1970]  has  derived  the  ap- 
propriate expressions  for  the  Rayleigh  waves  generated  by 
simple  sources  at  depth  in  a multilayered  halfspace. 
Harkrider,  et  al.  [1974]  extended  this  theory  to  the  case 
of  an  overpressure  applied  at  the  free  surface.  The  surface 
waves  generated  by  the  free  surface  boundary  condition  of 
interest  here  can  be  obtained  from  appropriate  modification 
of  these  results. 

For  far-field  body  waves  we  could  follow  the  procedure 
used  by  Fuchs  [1966]  , Hudson  [1969]  and  Bache  and  Harkrider 
[1976]  for  multilayered  media,  modifying  the  boundary  condi- 
tions appropriately.  However,  it  will  be  sufficient  for  our 
purposes  to  compute  the  far-field  body  waves  in  a homogeneous 
half space.  Treatment  of  this  elementary  problem  is  as  in 

Ewing,  Jardetsky  and  Press  [1957,  Chapter  2]. 

Let  us  assume  that  the  surface  loading  consists  of 
an  axisymmetric  distribution  of  normal  and  shear  stress 
components.  The  geometry  is  shown  in  Figure  1.  If  we 
Fourier  transform  all  transient  quantities,  the  free  sur- 
face boundary  conditions  are: 


Pzz  (r, 0,O,oi) 


P (r,ui)  , 

0 


Pzr {r'6'0'w> 


Q (r,u>)  , 
0 


(1) 
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Figure  1.  The  geometry  and  coordinate  system  for  an  axisym- 
metric  distribution  of  normal  and  shear  stress 
applied  at  the  surface  of  a multilayered  halfspace. 


126 


R-3119 


where  P and  P are  the  Fourier  transformed  normal  and 
zz  zr 

shear  stresses  and  P , Q represent  applied  stress  distri- 

0 0 

butions . 

To  solve  the  problem  we  first  compute  the  seismic 
waves  for  impulse  loading.  The  result  is  the  Green's 
function  solution  and  the  total  response  may  then  be  ob- 
tained by  area  integration. 

E. 2 GREEN'S  FUNCTION  FORMULATION 

For  polar  coordinates  the  Green's  function  for  inhomo- 
geneous boundary  conditions  is  the  response  to  boundary 
loading  of  the  form  <5(r-r")/r  [Morse  and  Feshbach  [1953], 
Chapter  7] . A convenient  representation  for  our  purpose  is 
(e.g.,  Watson  [1966],  Section  14.3), 


° (r~-r-  - = / Jv(kr)  Jv(kr")  k dk  . (2) 

0 


For  the  normal  traction  boundary  condition  the  appropriate 
form  of  (2)  has  v = 0.  For  the  shear  traction  the  v = 1 
case  is  appropriate.  Therefore,  the  Green's  function  bound- 
ary conditions  are: 


P (r,w) 
0 


Q (r,w) 
o 


J (kr)  J (kr")  k dk 
o o 


J (kr)  J (kr ")  k dk 
i l 


(3) 
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The  solution  is  obtained  by  solving  for  the  seismic 
waves  due  to  the  normal  and  shear  traction  boundary  condi- 
tions (3)  applied  separately.  The  result  is  the  Green's 
functions,  GN(r/r")  and  GT(r/r'),  for  the  normal  and  shear 
tractions.  The  total  solution,  UT(co),  is  then  given  by  the 
spatial  integration: 


UT(w) 


P (r, w) G (r/r " ) r "dr ' 
a N 


Q (r,u)G  (r/r')r'dr'. 
o T 

(4) 


E.3  SURFACE  WAVES  - THE  GREEN'S  FUNCTION 


The  formalism  of  Harkrider  [1964]  is  used  to  compute 
the  surface  waves.  Since  the  problem  is  axisymmetric , the 
6 dependent  quantities  vanish.  The  important  steps  in  the 
analysis  are  summarized  below.  The  notations  used  are 
essentially  the  same  of  those  of  Harkrider  [1964]. 


Assume  a cylindrical  coordinate  system  with  origin 
at  the  free  surface.  The  geometry  is  shown  in  Figure  1. 
The  nonvanishing  displacements  and  stresses  in  the  layers 
are  related  to  potentials  by 


q1 (r, z) 


a^1 

3r 


a2^1 

3rTz  ' 


w1  (r,z)  = + k24'i  , 


(5) 


Pzz(r'z)  = 2^i 


ill!  + ilii  + k2  iil 

3z2  3 z ! 6i  32 


- X . k2  01  , 
l a . ' 

1 


Prz(r'2)  “ Wi 


2d2<pi  . 233<P1  2 dip1 

3z3r  az2ar  0i  3r 
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Define 


oo 

I1  (r,z)  = / qi 

•'a 


( r , z ; k ) dk  , 


00 

wi(r,z)  = y wi(r,z; 


k)  dk  , 


?k)  dk  , 


j1  (r,  z)  -J  4>i  (r , z ; 


oo 

ji(rfz)  = y ^(r^k)  dk  , 

0 

OO 

’zz(r'z)  = f pzz.  (r'z;k)  dk  ' 
J0  1 

00 

’rz(r,z)  ~ f Prz  (r'z;k)  dk  • 

j zi 


(6) 


Since  the  location  of  material  boundaries  depends  only 
on  r,  the  r and  z dependence  of  the  potentials  can  be  sep- 
arated: 


4»i  (r , z ; k) 


4>.  (z)  J (kr)  , 
i o 


^ (r , z ; k) 


iK  (z)  J (kr)  . 
x o 


(7) 
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(8) 


If  (5)  and  (6)  are  substituted  into  (4)  we  obtain: 

,»R. 

qi(r,z;k)  H (kr)  , 

. V 

"i'r'2!k>  I - i J (kr)  , 


P (r , z ; k)  5 a (z)  J (kr)  , 

Z Z . K . q 

1 1 0 


P (r , z ; k)  = -i  t (z)  J (kr)  , 
rzi  Ri  i 


where  the  newly  introduced  normalized  velocities  and  stresses 
are  defined  in  terms  of  the  potentials  in  Eq.  (7)  of  Harkrider 
[1964]  . 

The  potentials  satisfy  the  wave  equations  and  there- 
fore may  be  written 


-ikrc  2 


ikr  z 

a . 


<J> . (z)  = A;  e 1 + A7"  e 1 


-ikrg  z ikr^  z 

<p^(z)  = u>'  e ^ + (1)7"  e ^ 


(9) 


Define 


-ikr  z . , 

1~1  - 


Air  " ~ka.  e 1 A£  ' 


2 ikraizi-l  _ 

Ar  - ~ka . e Ar  ' 

i 

ke.  ~ikr8 • zi-l 
= ik  -ji  e 1 Sj  , 


(10) 


A k8i  lkr0.zi-l 

uk'  = ik  -ji-  e 1 . 


130 


R-3119 


Then  continuity  of  stresses  and  displacements  at  layer  inter- 
faces may  be  expressed  by 


with  aR  a matrix  with  elements  given  by  Eq.  (16)  of 
Harkride^  [1964] . 

The  relationship  between  the  coefficients  in  (8)  and 
the  velocity-stress  vector  can  be  written  as 


R-3119 


. 

I 


I 


I 

l 

{ 


Applying  the  boundary  conditions  at  the  free  surface. 


• 

• 

UR 

UR 

— - (z  . ) 
c n-1 

2. 

c 

• 

• 

w„ 

OJ  _ 

R 

-r-  (z  ) 

c n-l 

R 

= A 

c 

°R  <zn-l> 

R 

a 

R 

n 

0 

TR  <*„-!> 

T 

R 

n 

0 

where  AR  = aR 


n-1 


n-2 


(13) 


From  (11)  and  (9)  together  with  the  condition  that  there 
be  no  radiation  from  infinite  depth. 


• 

A 

A " 
n 

UR 

FT  <*„-!> 

• 

l' 

WR 

FT  <zn-l> 

n 

U) 

n 

n 

aRn  <*n-l} 
n 

A 

0) 

n 

tR  <zn-l> 
n 
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Therefore,  we  get  Harkrider's  [1964]  equation  (62): 


— — 

— — 

A . 

A' 

n 

A _ 

w 

A 

n 

A - 

= j 

X 

CO 

n 

A 

Y 

* 

U) 

n 

Z 

— — 

_ 

with  J = E_  A_  and  W,  x,  Y,  Z the  elements  of  the  vector  on 
Kn  R 

the  right  side  of  (12);  i.e.,  the  free  surface  velocities  and 
stresses . 

The  system  of  equations  (14)  can  be  solved  for  W and 
X.  The  result  is 


(GN-LH)Y  + 


F 


R 


(RN-SL) Z 

t 


(RN-SL) Y + (RN-SK) Z 
F„ 


(16) 


where  FR,  G,  N,  etc.  are  combinations  of  elements  of  the  J 
matrix  given  by  Harkrider.  Rewrite  these  as 


(17) 
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By  (5)  and  (7) 


r WR 

wo  = u'  (r,0)  = - j ±-± 


J (kr)  dk  , 
0 


oo  U 


qo  = q 


'<r,0,  = -f  i 


— *•  J (kr)  dk 
c l 


Then,  using  (16), 


00 

/N 

* 


J (kr)  dk  , 
kFR  0 


o°  ( 3 ) 


q0  = -/  iEFT-  Vkr)  dk  • 

0 R 


For  the  surface  waves  we  evaluate  the  residue  contribution 
to  these  integrals.  For  ^ fixed  we  find  for  the  mode 

that  the  displacements  are: 


{w  }r  = j—  H(2)(k  r)  , 

3 ^ 3 


o R.  kn  75F 


j R 


‘ 


H(2)  (kR  r)  , 


where  (9FD/3k) .,  N^1^,  N*3'  are  evaluated  at  (w,kn  ) such 

K ] Kj  Kj  Kj 

that  F (w,kD  ) = 0.  These  can  be  combined  to  give 

R K-s 
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{q.}R  = - i 

o Rj 


u0 

w 


{ W } 


Hx(2)  (kR.r) 


H . 
3 


o R.  H ( 2 ) (k  r) 
o Rj 


(21) 


where  [u  /w  ]„  = -K/L,  the  ellipticity. 

o o « j 

It  remains  to  cast  (19)  in  a form  suitable  for  com- 


putation. We  have 


-(GN  - LN)  (Y  + 5.  Z)  , 


= -i (G*N  - L*N) 


where  G = iG*,  etc.  Also, 


(22) 


(23) 


Then  the  vertical  and  radial  components  of  the  Rayleigh  wave 
displacement  are 


{ 


(24) 


{qo} 

0 R. 


3 


kR. 

3 


- • 

u 

_! 

A_ 

V- 

- • 

U 

L 

.} 

• 

w 

0 

~R . 
3 

Hj 

. 

w 

L 0 J 

) 

H 


(2) 


<kR.r) 

3 
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where 


(G*N  - L*H) 


To  complete  the  Green's  function  formulation  it  is  neces- 
sary to  interpret  Y and  Z in  terms  of  the  boundary  con- 
ditions given  by  (3) . From  (5)  and  (7) 


[P  ] 
1 zz 


00 

-/ 

r=0 


J (kr) 
o 


dk 


fprz] 


z=0 


J (kr)  dk  , 
i 


(25) 


where  aR  = Y and  tr  = Z. 
o o 

Comparing  (24)  with  (1)  and  (3),  we  see  that 

a = k J (kr")  , 

Ro 

(26) 

tr  = ik  J (kr")  , 

R 1 

0 

give  the  desired  Green's  functions.  For  the  two  cases  of 
normal  and  shear  traction  they  are: 

Case  1 Applied  normal  traction:  The  Green's  function  is 

Gf.(r|r")  = -iir  An  J (kB  r")  H(2)(k_  r)  (27) 

N o 0 Rj 
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and  the  surface  wave  displacements  are  given  by 


(w  } = G* (r| r')  , 

8 R.  N 

D 

(28) 

N Hi(2)  (kR  r) 

{q  } = e.  -Ur — -3 — G®(r|r-)  , 

0 R.  3 H 2 (k_  r) 

j 0 Kj 

where  e.  = -i [u  /w  ]„  has  been  introduced. 

3 o o Hj 

Case  2 Applied  shear  traction:  The  Green's  function  is 

G®(r|r')  = -iir  AR  Ej  (kR  r')  H^2)  (kR  r)  , (29) 

and  the  displacements  are  given  by  (27)  with  Gf,(r|r')  re- 

S | „ 1 

placing  G^ ( r 1 r ) . 

E . 4 SURFACE  WAVES  - SOLUTION  FOR  A DISTRIBUTED  SURFACE 
LOAD 

Equations  (26)  and  (28)  give  the  Green's  functions  for 
normal  and  shear  traction  loading  applied  at  the  point  r = r'. 
From  these  Green's  functions  we  can  derive  the  solutions  for 
loads  distributed  over  the  free  surface.  The  basic  equation 
relating  the  Green's  function  to  the  total  solution  is  (4). 

Case  1:  Applied  Normal  Traction 


We  consider  the  boundary  conditions  (1)  and  assume 

initially  that  Q (r,uj)  vanishes.  Let  us  also  assume  that 

o 

the  loading  is  restricted  in  spatial  extent;  that  is, 

P (r,0,O,w)  = P (r,w)  , a < r < b, 

ZZ  Q “ — 


= 0 , elswhere. 
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Then 


and 


/ \N 
(w  } 

o Rj 


= f Po(r,o))  GN(r|r')  r'dr', 


(31) 


= -i-rrA  H(2^  (k  r)  f P (r',a))J  (k  r')r'dr't 
~Kj  o Kj  J o o Kj 


H!(2)  (kR.r) 

{q-}R,  = ei  Z(2)  j {W  } 


H (kD  r) 


o N 


Case  2:  Applied  Shear  Traction 

In  this  case  we  compute  the  displacements  for  boundary 

conditions  (1)  with  P (r,ai)  vanishing.  Again,  let  us  assume 

o 

that  the  loading  is  applied  for  a < r < b.  Then 


(w  = f Q (r ,o>)  G ( r | r ' ) r'dr' 
o Rj  J o T 


(32) 


* -i’?Rj  Ej  k,<2>  <kR 


r)  f Q (r',w)  J (k_  r')  r'dr', 
j 7 o ‘ Rj 


Hi(2)(kRr) 

{qa}R  = ei  “T75 3 — {<Vs  * 

0 Rj  3 HU)  (k  r)  0 s 

i ^ 


From  (31)  and  (32)  we  see  that  the  boundary  stress-time  his- 
tories appear  only  in  the  spatial  integrals 
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P(w)  = J P (r',to)  J (k  r')  r'dr' , 
a 0 o Rj 


Q(oj)  = / Q (r',uj)  J (k  r ) r'dr', 

•/  0 i K . 

a : 


(33) 


It  is  necessary  to  evaluate  these  integrals  for  each  fre- 
quency, w,  of  interest  for  the  surface  wave  calculations. 

In  general,  this  is  done  numerically. 

Since  the  problem  is  linear,  we  can  superimpose  the 
solutions  (31)  and  (32)  to  obtain  the  solution  when  both 
normal  and  shear  tractions  are  applied.  The  total  Rayleigh 
wave  displacement  is 


{W  } = {w  + (w  , 

o R • o R • o R • 

3 3 3 


«.v  - <Vr.  ♦ <vsRj 


(34) 


Substituting  from  (31)  and  (32)  , this  may  be  written 


{w  } = -iirA  H ( 2 } (k  r)  [P  + e . Q]  , 

0 Kj  ~Rj  0 Rj  3 


H(2)  (k  r) 

1 **2 

{qn}R  = ei  ~TT\ 3 — {w  }R  • 

0 Rj  3 H(2)  (k_  r)  0 Rj 

Rj 


(35) 


where  P and  Q are  given  by  the  integrals  (33)  . Finally,  note 
that  in  the  far-field  (kR  r >>  1)  we  can  use  the  asymptotic 
approximation  for  the  Hansel  functions  which  is 


Hf2)(Z)  * jn 

1 1 TTZ 


-i  (as  - - I) 

. 1 2 4 


(36) 


I 
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Then  the  expression  for  (q  } reduces  to 


{q  =ie.{w)_ 
o Rj  Do  Rj 


(37) 


For  the  fundamental  mode  (j  =0)  Rayleigh  wave  < 0 and 
(37)  prescribes  retrograde  elliptic  particle  motion. 


E. 5 BODY  WAVES  — THE  GREEN'S  FUNCTION 

As  we  pointed  out  in  the  introduction,  we  could  derive 
the  expressions  for  the  body  waves  due  to  distributed  surface 
loading  using  procedures  analogous  to  those  followed  for  sur- 
face waves.  However,  it  will  be  sufficient  to  consider  only 
the  homogeneous  halfspace  case  which  is  quite  easy  since  the 
solution  follows  directly  from  results  given  by  Ewing, 
Jardetsky  and  Press  [1957],  subsequently  referred  to  as  EJP. 

Case  1:  Applied  Normal  Traction 

Consider  the  Green's  function  boundary  conditions  (3) 
and  set  Q (r,u>)  = 0.  Then  from  EJP,  equations  (2-59)  to 
(2-68)  , it  is  quite  easy  to  show  that  the  diltational  poten- 
tial, may  be  written  as 


0 


N 


(2k2  - k2) 
F (k) 


-vz 


J (kr)J  (kr')  kdk, 
o o 


(38) 


where  F (k) , the  Rayleigh  function,  and  the  other  notations 
are  as  in  EJP. 

The  k-integral  (38)  represents  the  total  solution. 

The  far-field  body  wave  portion  of  the  solution  can  be 
shown  to  result  from  evaluating  (38)  in  the  complex  plane 
and  approximating  the  branch  line  contribution  to  the  integral 
by  using  the  method  of  steepest  descent.  A number  of  authors 
have  dealt  with  integrals  of  this  form  including  EJP,  Chapter 
2, and  Fuchs  [1966].  Bache  and  Harkrider  [1976,  Eq.  30]  write 
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the  solution  for  far-field  body  waves  in  compact  form.  Ob- 
serving this  formula,  the  solution  is 


i k (2k2  - k2) 

U F0T5  ra 


J (kr") 


-ik  R 
e a 
R 


(39) 


r 


2 

a 


c2/a2 


- 1, 


R2  = r2  + z2 


sine  = a/c  = r/R, 


where  6 is  the  takeoff  angle  measured  from  the  downward  verti- 
cal and  c is  horizontal  phase  velocity.  The  far-field  dis- 
placements are  then  found  by  differentiating  (39)  with  res- 
pect to  R and  retaining  only  terms  of  order  R-1.  The  result 
is  the  Green's  function 

(2k2  - k2)  “i^R 

(rlr'>  - ksk  — HIT)—  ra  J,  ,kr'>  I • (40) 


Case  2:  Applied  Shear  Traction 

Now  consider  (32)  and  set  P (r,w)  = o.  Then  following 
the  same  procedure  as  was  used  to  derive  (38)  , we  find  that 


♦s  - 7 / iw  e'vz  J0  <kr>  J,  <kr'>  dk-  <41> 

0 


The  Green's  function  for  far-field  body  waves  is  then 


Gg  (r  | r') 


2k  k 

- (V  2 

uF(k) 


-ik  R 

k8)  ra  (kr'}  I 


(42) 
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The  far-field  body  wave  Green's  functions  (40,  42)  are 
analogous  to  (27,  29)  for  surface  waves.  The  relation  be- 
tween the  Green's  functions  and  the  total  solution  for  the 
distributed  loads  (1)  is  given  by  (4) . Assuming  the  boundary 
loading  to  be  applied  in  the  region  a <_  r <_  b,  we  have  the 
total  far-field  body  wave  displacement 


h k (2k2  - k2)  'ikaR 

U(o)  =/  P (r',uO  k r J (kr')  | 

J o a UF  (k)  a „ R 


h k (k2-k2 ) e~lkatR 

+ f Vr''“>  2ka  ~ ikT~  rc  J,(kr')  R 


r 'dr  ' . 


This  may  be  reduced  to 

, , -ik  R . 

kk  r a - ~ j 

U(u>)  = uF^k)  ^ I | (2k2-k2 ) P(u>)  + (k2-k2)  Q ( ol»  ) | (44) 


where  P (w)  , Q(cj)  are  the  spatial  integrals  given  by  (33). 
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